Phương pháp ước lượng tối đa khả dĩ mục tiêu thường được gọi là TMLE là một công cụ mạnh mẽ trong phân tích nhân quả. Nhiều nhà nghiên cứu từng nghe về đặc tính song trùng bền vững của phương pháp này nhưng chỉ thực sự thấu hiểu nó khi tiến hành giả lập dữ liệu thực tế. Phương pháp này hoạt động cực kỳ hiệu quả khi một trong hai mô hình kết quả hoặc mô hình điều trị được thiết lập chính xác. Việc kết hợp thuật toán XGBoost cùng TMLE giúp tự động bắt trọn các mối quan hệ phức tạp trong dữ liệu mà không cần phải khai báo các tương tác thủ công. Bài viết này sẽ đi sâu vào cơ chế vận hành của phương pháp thông qua việc giả lập dữ liệu cụ thể trong môi trường R.
Khái Niệm Về Tmle
TMLE là một phương pháp thống kê tiên tiến được sử dụng để ước lượng các tác động nhân quả trong các nghiên cứu quan sát và thử nghiệm lâm sàng. Phương pháp này kết hợp linh hoạt giữa các thuật toán học máy và kỹ thuật thống kê truyền thống nhằm mang lại các ước lượng vững cho hiệu quả tác động của can thiệp, đồng thời kiểm soát tốt các yếu tố nhiễu.
Quy trình vận hành của TMLE gồm hai giai đoạn chính. Đầu tiên, hệ thống sẽ ước lượng mô hình kết quả và mô hình điều trị. Sau đó, các ước lượng này được sử dụng để hiệu chỉnh nhằm hướng trực tiếp đến tham số mục tiêu cần nghiên cứu. Cách tiếp cận này đặc biệt hữu ích trong các bối cảnh mà phương pháp truyền thống dễ bị lệch hoặc hoạt động kém hiệu quả do sự xuất hiện của các mối quan hệ phi tuyến phức tạp.
Ý Nghĩa Của Đặc Tính Song Trùng Bền Vững
Một ước lượng được coi là song trùng bền vững khi nó duy trì tính nhất quán nếu mô hình gán điều trị hay còn gọi là điểm xu hướng hoặc mô hình kết quả được xác định chính xác. Nói cách khác, bạn có hai cơ hội để đạt được một ước lượng nhân quả hợp lệ. Ngay cả khi một trong hai mô hình bị thiết lập sai, chỉ cần mô hình còn lại đúng, kết quả ước lượng cuối cùng vẫn đảm bảo độ tin cậy. Đặc tính này mang lại lợi thế vượt trội trong các nghiên cứu quan sát thực tế, nơi ranh giới của việc thiết lập đúng mô hình thường rất mơ hồ. Để đánh giá thuộc tính này, hai chỉ số quan trọng cần được xem xét là độ lệch bias và phương sai variance.
Phân Biệt Độ Lệch Bias Và Phương Sai Variance
Độ lệch bias và phương sai variance là hai nguồn sai số cơ bản trong mọi mô hình dự báo và thống kê. Độ lệch đại diện cho sai số hệ thống xảy ra khi mô hình liên tục đánh giá cao hơn hoặc thấp hơn giá trị thực tế của tham số. Độ lệch cao dễ dẫn đến hiện tượng học chưa đủ, khiến mô hình bỏ lỡ các quy luật ẩn sâu trong dữ liệu.
Ngược lại, phương sai thể hiện sự dao động của các dự báo trên các tập dữ liệu huấn luyện khác nhau. Phương sai cao dẫn đến hiện tượng quá khớp, nơi mô hình học cả những nhiễu ngẫu nhiên thay vì tín hiệu thực tế.
Công thức tính độ lệch bias được biểu diễn như sau:
Bias của theta_hat = E[theta_hat] - theta
Trong đó theta_hat là ước lượng của tham số theta và E[theta_hat] là giá trị kỳ vọng của ước lượng.
Đoạn mã giả lập bằng R dưới đây minh họa cách tính toán độ lệch:
1predicted_theta <- vector(mode = "numeric", length=1000)
2for (i in 1:1000) {
3 training_data <- dplyr::slice_sample(original_training_data, n = nrow(original_training_data), replace = T)
4 model <- glm(outcome~treatment+confounder,data=training_data)
5 outcome_hat_1 <- predict(model,newdata = training_data |> mutate(treatment = 1))
6 outcome_hat_0 <- predict(model,newdata = training_data |> mutate(treatment = 0))
7 predicted_theta[i] <- mean(outcome_hat_1) - mean(outcome_hat_0)
8}
9bias <- mean(predicted_theta) - thetaCông thức tính phương sai variance được biểu diễn như sau:
Var của theta_hat = E[(theta_hat - E[theta_hat])^2]
Đoạn mã giả lập bằng R tương ứng để tính toán phương sai:
1predicted_theta <- vector(mode = "numeric", length=1000)
2for (i in 1:1000) {
3 training_data <- dplyr::slice_sample(original_training_data, n = nrow(original_training_data), replace = T)
4 model <- glm(outcome~treatment+confounder,data=training_data)
5 outcome_hat_1 <- predict(model,newdata = training_data |> mutate(treatment = 1))
6 outcome_hat_0 <- predict(model,newdata = training_data |> mutate(treatment = 0))
7 predicted_theta[i] <- mean(outcome_hat_1) - mean(outcome_hat_0)
8}
9variance <- mean((predicted_theta-mean(predicted_theta))^2)
10# Hoặc sử dụng hàm có sẵn
11variance <- var(predicted_theta)Giả Lập Dữ Liệu Thực Tế
Hãy tiến hành xây dựng một tập dữ liệu giả lập có chủ đích với các mối quan hệ phi tuyến và tương tác phức tạp để kiểm tra tính năng của các phương pháp.
1library(tidyverse)
2set.seed(1)
3n <- 10000
4W1 <- rnorm(n)
5W2 <- rnorm(n)
6W3 <- rbinom(n, 1, 0.5)
7W4 <- rnorm(n)
8# Mô hình điểm xu hướng thực tế
9A <- rbinom(n, 1, plogis(-0.5 + 0.8*W1 + 0.5*W2^2 + 0.3*W3 - 0.4*W1*W2 + 0.2*W4))
10# Mô hình kết quả thực tế
11Y <- rbinom(n, 1, plogis(-1 + 0.2*A + 0.6*W1 - 0.4*W2^2 + 0.5*W3 + 0.3*W1*W3 + 0.2*W4^2))
12# Tính toán giá trị ATE thực tế
13logit_Y1 <- -1 + 0.2 + 0.6*W1 - 0.4*W2^2 + 0.5*W3 + 0.3*W1*W3 + 0.2*W4^2
14logit_Y0 <- -1 + 0 + 0.6*W1 - 0.4*W2^2 + 0.5*W3 + 0.3*W1*W3 + 0.2*W4^2
15Y1_true <- plogis(logit_Y1)
16Y0_true <- plogis(logit_Y0)
17true_ATE <- mean(Y1_true - Y0_true)
18df <- tibble(W1 = W1, W2 = W2, W3 = W3, W4 = W4, A = A, Y = Y,
19 true_ATE = true_ATE, Y1_true = Y1_true, Y0_true = Y0_true)Giá trị hiệu quả tác động trung bình thực tế thu được từ dữ liệu này là 0.0373518. Chúng ta sẽ dùng giá trị mốc này để đối chiếu hiệu năng của các mô hình kế tiếp.
Xây Dựng Các Mô Hình Ước Lượng Từ Sai Đến Đúng
Mô hình kết quả sai lệch
Nếu chỉ sử dụng một mô hình hồi quy tuyến tính thông thường không khai báo các tương tác bậc hai:
1model <- glm(Y ~ A + W1 + W2 + W3 + W4, family = "binomial")
2summary(model)Hệ số của biến điều trị A ước lượng được là âm 0.0501418, hoàn toàn ngược lại với tác động thực tế của dữ liệu.
Hàm g-computation ứng dụng
Xây dựng hàm tính toán hiệu quả tác động trung bình dựa trên mô hình đã khớp:
1g_comp <- function(model,data,ml=F) {
2 if (ml==T) {
3 y1 <- predict(model, new_data=data |> mutate(A=as.factor(1)), type = "prob")[,2] |> pull()
4 y0 <- predict(model, new_data=data |> mutate(A=as.factor(0)), type = "prob")[,2] |> pull()
5 } else {
6 y1 <- predict(model, newdata=data |> mutate(A=1), type = "response")
7 y0 <- predict(model, newdata=data |> mutate(A=0), type = "response")
8 }
9 return(mean(y1-y0))
10}
11g_comp(model,df)Kết quả trả về là âm 0.009823307, xác nhận sự sai lệch nghiêm trọng khi bỏ qua các tương tác phi tuyến.
Mô hình kết quả chính xác
Khi khai báo đầy đủ các tương tác và biến số bậc hai thực tế vào mô hình:
1model <- glm(Y ~ A + W1 + I(W2^2) + W3 + W1:W2 + W4, family = "binomial")
2g_comp(model,df)Kết quả trả về đạt mức 0.03576854, rất sát với giá trị thực tế của ATE.
Mô hình điều trị sai lệch dùng phương pháp nghịch đảo xác suất gán điều trị
Sử dụng phương pháp nghịch đảo xác suất gán điều trị với mô hình gán điều trị bị thiết lập sai:
1ps_model <- glm(A ~ W1 + W2 + W3 + W4, family = "binomial")
2ps <- ps_model$fitted.values
3ps_final <- pmax(pmin(ps, 0.95), 0.05)
4weights <- ifelse(A == 1, 1/ps_final, 1/(1-ps_final))
5model <- glm(Y ~ A, family = "binomial", weights = weights)
6g_comp(model,df)Kết quả thu được tiếp tục sai lệch nặng nề với giá trị âm 0.0095777.
Mô hình điều trị chính xác
Khi mô hình gán điều trị được khai báo đúng các thành phần phi tuyến:
1ps_model <- glm(A ~ W1 + I(W2^2) + W3 + W1:W2 + W4, family = "binomial")
2ps <- ps_model$fitted.values
3ps_final <- pmax(pmin(ps, 0.95), 0.05)
4weights <- ifelse(A == 1, 1/ps_final, 1/(1-ps_final))
5model <- glm(Y ~ A, family = "binomial", weights = weights)
6g_comp(model,df)Kết quả cải thiện rõ rệt đạt mức 0.03874379, tiệm cận rất gần với giá trị mục tiêu.
Giải pháp học máy với XGBoost
Trong thực tế, chúng ta hầu như không thể biết trước cấu trúc phi tuyến phức tạp của dữ liệu để khai báo thủ công. Sử dụng thuật toán XGBoost để tự động dò tìm cấu trúc này:
1library(tidymodels)
2library(doParallel)
3library(future)
4workers <- parallel::detectCores(logical = FALSE) - 1
5plan(multisession, workers = workers)
6all_cores <- parallel::detectCores(logical = FALSE)
7cl <- makePSOCKcluster(all_cores - 1)
8registerDoParallel(cl)
9df_ml <- df |>
10 select(Y,A,W1,W2,W3,W4) |>
11 mutate(Y = as.factor(Y),
12 A = as.factor(A))
13xgb_spec <- boost_tree(
14 trees = 1000,
15 tree_depth = tune(),
16 min_n = tune(),
17 loss_reduction = tune(),
18 mtry = tune(),
19 learn_rate = tune()
20) %>%
21 set_engine("xgboost") %>%
22 set_mode("classification")
23xgb_wf <- workflow() %>%
24 add_model(xgb_spec) %>%
25 add_formula(Y ~ .)
26xgb_grid <- grid_space_filling(
27 tree_depth(),
28 min_n(),
29 loss_reduction(),
30 finalize(mtry(),df_ml),
31 learn_rate(),
32 size = 20
33)
34set.seed(1)
35folds <- vfold_cv(df_ml, v = 5)
36xgb_res <- tune_grid(
37 xgb_wf,
38 resamples = folds,
39 grid = xgb_grid,
40 control = control_grid(save_pred = TRUE,
41 parallel_over = "everything")
42)
43best_xgb <- select_best(xgb_res, metric = "roc_auc")
44final_xgb <- finalize_workflow(xgb_wf, best_xgb)
45final_fit <- fit(final_xgb, data = df_ml)
46g_comp(final_fit, df_ml, T)Kết quả trả về là 0.03447109. XGBoost đã tự động xử lý các mối quan hệ phức tạp và cho ra ước lượng rất chính xác mà không cần bất kỳ sự can thiệp thủ công nào vào công thức thiết lập.
Quy Trình Thực Hiện TMLE Chi Tiết
Hãy cùng xem xét quy trình TMLE từng bước một trong trường hợp mô hình kết quả bị thiết lập sai nhưng mô hình điều trị được thiết lập đúng.
Bước 1: Khởi tạo mô hình kết quả
1model_outcome <- glm(Y ~ A + W1 + W2 + W3 + W4, family = "binomial") # Thiết lập sai
2model_outcome_all <- predict(model_outcome, newdata = df, type = "response")
3model_outcome_1 <- predict(model_outcome, newdata = df |> mutate(A = 1), type = "response")
4model_outcome_0 <- predict(model_outcome, newdata = df |> mutate(A = 0), type = "response")Bước 2: Thiết lập mô hình điều trị và tính biến đồng biến thông minh
1model_treatment <- glm(A ~ W1 + I(W2^2) + W3 + W1:W2 + W4, family = "binomial") # Thiết lập đúng
2ps <- model_treatment$fitted.values
3ps_final <- ps
4a_1 <- 1/(predict(model_treatment, df, type = "response"))
5a_0 <- -1/(1 - predict(model_treatment, df, type = "response"))
6clever_covariate <- ifelse(A == 1, 1/ps_final, -1/(1-ps_final))Bước 3: Ước lượng tham số biến động epsilon
1epsilon_model <- glm(Y ~ -1 + offset(qlogis(model_outcome$fitted.values)) + clever_covariate, family = "binomial")
2epsilon <- epsilon_model$coefficientsBước 4: Cập nhật kết quả ban đầu dựa trên thông tin hiệu chỉnh
1updated_outcome_1 <- plogis(qlogis(model_outcome_1)+epsilon*a_1)
2updated_outcome_0 <- plogis(qlogis(outcome_0 <- plogis(qlogis(model_outcome_0)+epsilon*a_0))Bước 5: Tính toán hiệu quả tác động trung bình ATE mới
1ate <- mean(updated_outcome_1-updated_outcome_0)
2print(ate)Kết quả đạt mức 0.04068321, kéo giá trị ước lượng từ mức sai lệch nghiêm trọng ban đầu của mô hình kết quả về sát với thực tế nhờ vào mô hình điều trị chính xác.
Bước 6: Ước lượng sai số chuẩn và khoảng tin cậy
1se <- sqrt(var((Y-model_outcome_all)*clever_covariate+updated_outcome_1-updated_outcome_0-ate)/n)
2pval <- 2*(1-pnorm(abs(ate/se)))
3print(paste0("ATE: ",round(ate,3), " [95%CI ", round(ate-1.96*se,3),"-",round(ate+1.96*se,3),", p=",round(pval,3),"]"))Đánh Giá Và So Sánh Toàn Diện Các Mô Hình
Để kiểm chứng tính ổn định, chúng ta tiến hành lấy mẫu lại 1000 lần với quy mô mẫu là 6000 quan sát để so sánh bốn hướng tiếp cận khác nhau.
Bảng so sánh hiệu năng chi tiết dựa trên 1000 lần chạy giả lập:
| Tên Mô Hình | Độ Lệch Bias | Phương Sai Variance |
| --- | --- | --- |
| Mô hình hồi quy kết quả đúng | -0.0009778 | 0.0001410 |
| Mô hình điều trị đúng với trọng số IPW | 0.0019444 | 0.0001581 |
| Mô hình hồi quy kết quả sai lệch | -0.0464683 | 0.0001354 |
| Mô hình thuần XGBoost cho kết quả | -0.0164066 | 0.0000459 |
| Mô hình TMLE kết hợp XGBoost cỡ mẫu tìm kiếm 5 | 0.0066270 | 0.0001399 |
| Mô hình TMLE kết hợp XGBoost cỡ mẫu tìm kiếm 20 | 0.0041615 | 0.0001645 |
Kết quả so sánh cho thấy khi các mô hình hồi quy truyền thống được xác định chính xác, chúng cho độ lệch thấp nhất. Tuy nhiên, trong điều kiện thực tế khó xác định đúng dạng hàm, mô hình TMLE kết hợp với việc tối ưu hóa siêu tham số của XGBoost mang lại sự cân bằng tối ưu khi giảm thiểu độ lệch xuống mức rất thấp, đồng thời kiểm soát phương sai ở mức chấp nhận được tương đương với hồi quy chuẩn.
Giải Pháp Thay Thế Trong Thống Kê Truyền Thống
Bên cạnh TMLE, mô hình ước lượng nghịch đảo xác suất gán điều trị tăng cường hay còn gọi là AIPW cũng là một đại diện tiêu biểu cho trường phái song trùng bền vững. AIPW kết hợp mô hình điểm xu hướng và mô hình kết quả để đưa ra ước lượng tối ưu và nhất quán chỉ cần một trong hai mô hình đúng. Trong thực tế, các nhà nghiên cứu có thể sử dụng các gói thư viện chuyên dụng của AIPW trong R để thay thế hoặc đối chiếu trực tiếp với kết quả của TMLE.
✨ Phương pháp TMLE chứng minh giá trị vượt trội khi kết hợp sức mạnh khai phá dữ liệu tự động của các thuật toán học máy hiện đại với nền tảng lý thuyết vững chắc của thống kê nhân quả, giúp các nhà khoa học dữ liệu tự tin đưa ra các kết luận chính xác hơn mà không lo ngại về sai lệch do thiết lập sai dạng hàm của mô hình.
Câu hỏi tư duy: Trong các nghiên cứu quan sát có số lượng biến nhiễu cực kỳ lớn và xuất hiện hiện tượng đa cộng tuyến mạnh, liệu việc sử dụng thuật toán học máy phức tạp cho cả mô hình điều trị lẫn mô hình kết quả trong TMLE có làm tăng nguy cơ quá khớp và ảnh hưởng thế nào đến phương sai của ước lượng ATE cuối cùng?
