library(dplyr)
set.seed(123)
n <- 100 # RCT 대상자 수
baseline_SBP <- rnorm(n, 150, 20) # 기저 혈압
treat_num <- rbinom(n, 1, 0.5) # 무작위배정
treat <- factor(treat_num, levels=c(0,1))
lp0 <- -2 + 0.1*(baseline_SBP-150)
lp1 <- lp0 - 1.5
p0 <- plogis(lp0)
p1 <- plogis(lp1)
y <- rbinom(n,1,ifelse(treat_num==1, p1, p0))
df <- data.frame(y, treat, baseline_SBP)
calc_measures <- function(d) {
# Unadjusted
r1_un <- mean(d$y[d$treat=="1"])
r0_un <- mean(d$y[d$treat=="0"])
rd_un <- r1_un - r0_un
# g-computation
m_gc <- glm(y ~ treat + baseline_SBP, family=binomial, data=d)
d1 <- transform(d, treat=factor(1,levels=c(0,1)))
d0 <- transform(d, treat=factor(0,levels=c(0,1)))
r1_gc <- mean(predict(m_gc,newdata=d1,type="response"))
r0_gc <- mean(predict(m_gc,newdata=d0,type="response"))
rd_gc <- r1_gc - r0_gc
rr_gc <- r1_gc / r0_gc
or_gc <- (r1_gc/(1-r1_gc)) / (r0_gc/(1-r0_gc))
# IPTW
ps_mod <- glm(as.numeric(as.character(treat)) ~ baseline_SBP,
family=binomial, data=d)
ps <- predict(ps_mod,type="response")
w <- ifelse(d$treat=="1", 1/ps, 1/(1-ps))
r1_ip <- sum(w * d$y * (d$treat=="1")) / sum(w * (d$treat=="1"))
r0_ip <- sum(w * d$y * (d$treat=="0")) / sum(w * (d$treat=="0"))
rd_ip <- r1_ip - r0_ip
rr_ip <- r1_ip / r0_ip
or_ip <- (r1_ip/(1-r1_ip)) / (r0_ip/(1-r0_ip))
c(rd_un, r1_gc, r1_ip,
rd_gc, rr_gc, or_gc,
rd_ip, rr_ip, or_ip)
}
# 3) 원본 지표 계산
orig <- calc_measures(df)
# 4) 부트스트랩으로 SE 계산 (B=200)
B <- 200
boot_vals <- replicate(B, {
idx <- sample(seq_len(n), n, replace=TRUE)
calc_measures(df[idx,])
})
se_vals <- apply(boot_vals, 1, sd)
# 5) 결과 정리
res <- data.frame(
Method = c("Unadjusted","G-computation","IPTW"),
RiskDiff = c(orig[1], orig[4], orig[7]),
SE_RiskDiff = c(se_vals[1],se_vals[4],se_vals[7]),
RR = c(NA, orig[5], orig[8]),
SE_RR = c(NA, se_vals[5],se_vals[8]),
OR = c(NA, orig[6], orig[9]),
SE_OR = c(NA, se_vals[6],se_vals[9])
)
print(res)
