必要な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 errorslibrary(readxl)
library(tableone)
library(lubridate)
library(skimr)
library(svglite)
library(haven)
library(summarytools)##
## 次のパッケージを付け加えます: 'summarytools'
##
## 以下のオブジェクトは 'package:tibble' からマスクされています:
##
## viewlibrary(naniar)##
## 次のパッケージを付け加えます: 'naniar'
##
## 以下のオブジェクトは 'package:skimr' からマスクされています:
##
## n_completelibrary(devtools)## 要求されたパッケージ usethis をロード中ですlibrary(reader)## 要求されたパッケージ NCmisc をロード中です
##
## 次のパッケージを付け加えます: 'reader'
##
## 以下のオブジェクトは 'package:NCmisc' からマスクされています:
##
## cat.path, get.ext, rmv.extlibrary(stringr)
library(missForest)
library(mice)##
## 次のパッケージを付け加えます: 'mice'
##
## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter
##
## 以下のオブジェクトは 'package:base' からマスクされています:
##
## cbind, rbindlibrary(geepack)
library(ggplot2)
library(cobalt)## cobalt (Version 4.4.1, Build Date: 2022-11-03)library(tableone)
library(sandwich)
library(survey)## 要求されたパッケージ grid をロード中です
## 要求されたパッケージ Matrix をロード中です
##
## 次のパッケージを付け加えます: 'Matrix'
##
## 以下のオブジェクトは 'package:tidyr' からマスクされています:
##
## expand, pack, unpack
##
## 要求されたパッケージ survival をロード中です
##
## 次のパッケージを付け加えます: 'survey'
##
## 以下のオブジェクトは 'package:graphics' からマスクされています:
##
## dotchartlibrary(WeightIt)
library(twang)## To reproduce results from prior versions of the twang package, please see the version="legacy" option described in the documentation.library(naniar)
pacman::p_load(norm2, miceadds, lmtest, car, ROCR, pROC,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)
# sATEのためのstabilized weightsを計算
treatment_prob <- mean(filtered_data$antibio_exposure)
stabilized_weights <- ifelse(filtered_data$antibio_exposure == 1, treatment_prob/ps, (1-treatment_prob)/(1-ps))
# Use the stabilized weights to fit the outcome model
iptw_model <- glm(outcome ~ antibio_exposure, data = filtered_data, weights = stabilized_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!# 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_sATE <- data.frame(
Child_Score = "a,b",
OR = or,
Lower_CI = lower,
Upper_CI = upper,
P_value = p
)
print(result_a_b_sATE)## Child_Score OR Lower_CI Upper_CI P_value
## antibio_exposure a,b 0.9009794 0.4168081 1.947572 0.7909115c
library(dplyr)
# `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)
# sATEのためのstabilized weightsを計算
treatment_prob_c <- mean(imputed_data_c$antibio_exposure)
stabilized_weights_c <- ifelse(imputed_data_c$antibio_exposure == 1, treatment_prob_c/ps_c, (1-treatment_prob_c)/(1-ps_c))
# Use the stabilized weights to fit the outcome model
iptw_model_c <- glm(outcome ~ antibio_exposure, data = imputed_data_c, weights = stabilized_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_sATE <- data.frame(
Child_Score = "c",
OR = or_c,
Lower_CI = lower_c,
Upper_CI = upper_c,
P_value = p_c
)
print(result_c_sATE)## Child_Score OR Lower_CI Upper_CI P_value
## antibio_exposure c 1.905181 0.9607718 3.777917 0.0649754library(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 2a ステロイド、nsaids、aptt cateは抜かないと回らない。
library(dplyr)
# `child_score` をフィルタリングして 'a' を取得
imputed_data_a <- imputed_data %>% filter(child_score == "a")
# Propensity score modelの適合
propensity_model_a <- 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_a,
family = binomial("logit"))
# Propensity scoresを計算
ps_a <- predict(propensity_model_a, type = "response", newdata = imputed_data_a)
# sATEのためのstabilized weightsを計算
treatment_prob_a <- mean(imputed_data_a$antibio_exposure)
stabilized_weights_a <- ifelse(imputed_data_a$antibio_exposure == 1, treatment_prob_a/ps_a, (1-treatment_prob_a)/(1-ps_a))
# Use the stabilized weights to fit the outcome model
iptw_model_a <- glm(outcome ~ antibio_exposure, data = imputed_data_a, weights = stabilized_weights_a, family = binomial("logit"))## Warning in eval(family$initialize): non-integer #successes in a binomial glm!# 結果を取得
estimate_a <- coef(iptw_model_a)[2]
std_error_a <- coef(summary(iptw_model_a))[2, "Std. Error"]
conf_int_a <- confint(iptw_model_a, 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_a <- exp(estimate_a)
lower_a <- exp(estimate_a - 1.96 * std_error_a)
upper_a <- exp(estimate_a + 1.96 * std_error_a)
# p-valueを計算
p_a <- 2 * (1 - pnorm(abs(estimate_a / std_error_a)))
# 結果をdata frameに保存
result_a_sATE <- data.frame(
Child_Score = "a",
OR = or_a,
Lower_CI = lower_a,
Upper_CI = upper_a,
P_value = p_a
)
print(result_a_sATE)## Child_Score OR Lower_CI Upper_CI P_value
## antibio_exposure a 0.8665251 0.1023304 7.337659 0.8954255b
library(dplyr)
# `child_score` をフィルタリングして 'b' を取得
imputed_data_b <- imputed_data %>% filter(child_score == "b")
# Propensity score modelの適合
propensity_model_b <- 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_b,
family = binomial("logit"))
# Propensity scoresを計算
ps_b <- predict(propensity_model_b, type = "response", newdata = imputed_data_b)
# sATEのためのstabilized weightsを計算
treatment_prob_b <- mean(imputed_data_b$antibio_exposure)
stabilized_weights_b <- ifelse(imputed_data_b$antibio_exposure == 1, treatment_prob_b/ps_b, (1-treatment_prob_b)/(1-ps_b))
# Use the stabilized weights to fit the outcome model
iptw_model_b <- glm(outcome ~ antibio_exposure, data = imputed_data_b, weights = stabilized_weights_b, family = binomial("logit"))## Warning in eval(family$initialize): non-integer #successes in a binomial glm!# 結果を取得
estimate_b <- coef(iptw_model_b)[2]
std_error_b <- coef(summary(iptw_model_b))[2, "Std. Error"]
conf_int_b <- confint(iptw_model_b, 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_b <- exp(estimate_b)
lower_b <- exp(estimate_b - 1.96 * std_error_b)
upper_b <- exp(estimate_b + 1.96 * std_error_b)
# p-valueを計算
p_b <- 2 * (1 - pnorm(abs(estimate_b / std_error_b)))
# 結果をdata frameに保存
result_b_sATE <- data.frame(
Child_Score = "b",
OR = or_b,
Lower_CI = lower_b,
Upper_CI = upper_b,
P_value = p_b
)
print(result_b_sATE)## Child_Score OR Lower_CI Upper_CI P_value
## antibio_exposure b 0.7932464 0.3266856 1.926133 0.6088369interaction 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
# sATEのためのstabilized weightsを計算
treatment_prob <- mean(imputed_data$antibio_exposure)
stabilized_weights <- ifelse(imputed_data$antibio_exposure == 1, treatment_prob/ps, (1-treatment_prob)/(1-ps))
# stabilized weightsを用いてlogistic regression modelを適合
iptw_model <- glm(outcome ~ antibio_exposure * child_grouped, data = imputed_data, weights = stabilized_weights, 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.52"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)
# sATEのためのstabilized weightsを計算
treatment_prob <- mean(imputed_data$antibio_exposure)
stabilized_weights_interaction <- ifelse(imputed_data$antibio_exposure == 1, treatment_prob/ps_interaction, (1-treatment_prob)/(1-ps_interaction))
# 交互作用項を持つモデルを適合
iptw_interaction_model <- glm(outcome ~ antibio_exposure * child_score,
data = imputed_data,
weights = stabilized_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 = stabilized_weights_interaction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2540 -0.4482 -0.3550 -0.2917 3.6350
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.1749 0.5105 -6.219 4.99e-10 ***
## antibio_exposure -0.2997 1.0748 -0.279 0.78034
## child_scoreb 0.4933 0.5617 0.878 0.37976
## child_scorec 1.7854 0.5489 3.253 0.00114 **
## antibio_exposure:child_scoreb 0.4435 1.1511 0.385 0.70002
## antibio_exposure:child_scorec 0.7035 1.1277 0.624 0.53273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 529.99 on 789 degrees of freedom
## Residual deviance: 487.17 on 784 degrees of freedom
## AIC: 507.06
##
## Number of Fisher Scoring iterations: 6# 交互作用項のp値を取得
interaction_p_value <- coef(summary(iptw_interaction_model))["antibio_exposure:child_scorec", "Pr(>|z|)"]
print(paste("交互作用項のp値: ", interaction_p_value))## [1] "交互作用項のp値: 0.532727195257351"