Tại đây, tôi sẽ trình bày một ví dụ mô phỏng nhỏ trong R. Đầu tiên, chúng ta sử dụng phương pháp hai biến X1 và X2 phân phối chuẩn và Y mô phỏng một tập dữ liệu lớn, trong đó Y tuân theo mô hình logic đã cho X1 và X2.

Đầu tiên, chúng ta mô phỏng một tập dữ liệu hoàn chỉnh rất lớn:

# Mô phỏng tập dữ liệu hoàn chỉnh
expit <- function(x) {
  exp(x) / (1 + exp(x))
}

n <- 100000
x <- mvrnorm(n, mu = c(0, 0), Sigma = matrix(c(1, 0.2, 0.2, 1), nrow = 2))
x1 <- x[, 1]
x2 <- x[, 2]
y <- 1 * (runif(n) < expit(-3 + log(3) * x1 + log(3) * x2))
mean(y)

Tiếp theo, chúng ta ước tính tỷ lệ rủi ro biên khi thay đổi X1 từ 1 thành 0:

# Ước tính tỷ lệ rủi ro biên cho x1 = 1 so với x1 = 0, chuẩn hóa theo dữ liệu hoàn chỉnh
# Sau này cho MI, chúng ta sẽ viết một hàm lấy tập dữ liệu và trả về ước tính này
tinh_tyle_rui_ro_biên <- function(tap_du_lieu) {
  mo_hinh <- glm(y ~ x1 + x2, family = "binomial", data = tap_du_lieu)
  # Dự đoán rủi ro khi x1 = 0
  rui_ro_0 <- expit(coef(mo_hinh)[1] + coef(mo_hinh)[2] * 0 + coef(mo_hinh)[3] * tap_du_lieu$x2)
  # Dự đoán rủi ro khi x1 = 1
  rui_ro_1 <- expit(coef(mo_hinh)[1] + coef(mo_hinh)[2] * 1 + coef(mo_hinh)[3] * tap_du_lieu$x2)
  # Ước tính tỷ lệ rủi ro biên
  return(mean(rui_ro_1) / mean(rui_ro_0))
}

tap_du_lieu_hoan_chinh <- data.frame(y = y, x1 = x1, x2 = x2)
tinh_tyle_rui_ro_biên(tap_du_lieu_hoan_chinh)

Tiếp theo, chúng ta sử dụng một cơ chế mà Sullivan và các đồng nghiệp đã xem xét, tạo ra một số giá trị bị thiếu trong Y và X2:

# Tạo dữ liệu bị thiếu theo Sullivan và các đồng nghiệp
z1 <- x1 / 0.2^0.5
r_y <- 1 * (runif(n) < expit(2.5 + 2 * z1))
r_x2 <- 1 * (runif(n) < expit(2.5 + 2 * z1))
tap_du_lieu_quan_sat <- tap_du_lieu_hoan_chinh
tap_du_lieu_quan_sat$y[r_y == 0] <- NA
tap_du_lieu_quan_sat$x2[r_x2 == 0] <- NA

Bây giờ chúng ta có thể ước tính các giá trị bị thiếu trong Y và X2. Chỉ định mô hình kết quả logic cho các kết quả bị thiếu và các giá trị biến协变量 bị thiếu từ mô hình điền tương thích với mô hình kết quả logic:

num_imputations <- 10
tap_du_lieu_duoc_dien <- amelia(tap_du_lieu_quan_sat, m = num_imputations, 
                              idvars = c("y", "x1", "x2"), 
                              noms = c("y", "x2"), 
                             ords = c(), 
                              logs = c("y", "x2"))

Cuối cùng, chúng ta có thể áp dụng hàm đã định nghĩa trước đó để ước tính tỷ lệ rủi ro biên cho mỗi tập dữ liệu được điền, và kết hợp chúng bằng quy tắc Rubin (tức là lấy trung bình của logarit tỷ lệ rủi ro):

est_log_rr <- array(0, dim = num_imputations)
for (i in 1:num_imputations) {
  est_log_rr[i] <- log(tinh_tyle_rui_ro_biên(tap_du_lieu_duoc_dien$imputations[[i]]))
}
# Ước tính kết hợp của log tỷ lệ rủi ro là
mean(est_log_rr)
# và ước tính của tỷ lệ rủi ro
exp(mean(est_log_rr))

Chúng ta nhận được một ước tính rất gần với ước tính dữ liệu hoàn chỉnh sau khi điền.