M04 · Mô phỏng vi mô (microsimulation)

Khi mỗi người bệnh mang một lịch sử riêng: mô hình chuyển trạng thái cấp cá thể trong R

Vì sao bạn nên quan tâm?

Ở phòng khám ngoại trú HIV, người bệnh ổn định lĩnh thuốc TLD 3 tháng một lần, và đến hẹn thì được chỉ định xét nghiệm tải lượng virus. Theo bảng giá BHYT của một bệnh viện hạng đặc biệt, một lần xét nghiệm là 979,700 VNĐ, bằng khoảng 4 ngày lương tối thiểu vùng I và đắt hơn cả một năm thuốc TLD theo giá tham chiếu quốc tế (819,016 VNĐ, giá tại cảng đi) (Bộ Y tế 2024; The Global Fund 2026; Chính phủ 2025). Cuối năm 2024, cả nước có 184,214 người đang điều trị ARV, phần lớn lĩnh thuốc qua BHYT, và khoảng 96% đã ức chế được virus (Báo Lào Cai 2025). Nếu mỗi người xét nghiệm một lần mỗi năm theo mức giá này, khoản chi là khoảng 180 tỷ VNĐ. Module này dựng mô hình từng người bệnh để trả lời một câu hỏi ngân sách có thật: xét nghiệm tải lượng bao lâu một lần là đáng tiền?

Tổng quan

Mô hình Markov đoàn hệ ở M02 theo dõi tỷ lệ đoàn hệ trong từng trạng thái và không nhớ ai đã đi qua con đường nào. Module này dựng mô phỏng vi mô (microsimulation): cho từng người đi qua các trạng thái, mỗi người mang đặc điểm và lịch sử riêng, rồi lấy trung bình trên rất nhiều người (Siebert và c.s. 2012; Krijkamp và c.s. 2018). Ta viết cùng một mô hình theo ba cách: vòng lặp theo từng người, vectơ hóa, và C++ qua Rcpp. Ba cách cho kết quả giống hệt nhau đến từng chữ số, và khi tắt các phần riêng của mô phỏng vi mô thì nó tái tạo được mô hình đoàn hệ. Sau đó ta học cách chọn số cá thể, kiểm soát sai số Monte Carlo và chạy phân tích độ nhạy xác suất (PSA) lồng hai vòng.

Mục tiêu học tập

Sau module này, người học có thể:

  1. Nhận ra khi nào mô hình đoàn hệ không đủ: tính không đồng nhất có quan hệ phi tuyến, phụ thuộc lịch sử, bùng nổ số trạng thái.
  2. Viết một mô phỏng vi mô hoàn chỉnh trong R cơ bản, kèm điểm kiểm tra và dùng chung số ngẫu nhiên giữa các chiến lược.
  3. Tăng tốc mô hình bằng vectơ hóa, tính toán song song và Rcpp, đồng thời chứng minh mọi phiên bản cho cùng kết quả.
  4. Chọn số cá thể dựa trên sai số Monte Carlo, tách phương sai PSA thành phần tham số và phần nhiễu, và báo cáo minh bạch.

Kiến thức cần có

  • M00 Bài 6: số ngẫu nhiên, set.seed(), sai số Monte Carlo.
  • M02: trạng thái, chu kỳ, nguy cơ cạnh tranh (Bài 2), vết đoàn hệ và hiệu chỉnh trong chu kỳ (Bài 3), trạng thái đường hầm (Bài 4).
  • M03: khái niệm PSA, CEAC và EVPI.

Ví dụ xuyên suốt

Người trưởng thành nhiễm HIV đang điều trị ARV (antiretroviral therapy, ART) tại một phòng khám ngoại trú, dùng phác đồ bậc một (first-line regimen) TLD (tenofovir, lamivudin, dolutegravir) và đã ức chế được virus (viral suppression). Khi ổn định, họ tái khám và lĩnh thuốc 3 tháng một lần (Bộ Y tế 2021; Trần Lam 2023). Mức tuân thủ điều trị (adherence) khác nhau giữa các cá thể. Người tuân thủ kém dễ bị bật lại tải lượng virus (viral load, VL). Khi virus không còn bị ức chế (thất bại virus học, virological failure), kháng thuốc với phác đồ bậc một có thể xuất hiện, nhưng ở thời dolutegravir điều này hiếm (Matsumoto và c.s. 2025; World Health Organization 2024). Theo hướng dẫn quốc gia, người bệnh được xét nghiệm tải lượng 12 tháng một lần; khi tải lượng cao, họ được tư vấn tăng cường tuân thủ (tối đa 3 buổi), xét nghiệm lại sau 3 tháng, và chuyển sang phác đồ bậc hai (second-line regimen) AZT + 3TC + LPV/r nếu vẫn thất bại (Bộ Y tế 2021). Ta so sánh chiến lược hiện hành (12 tháng/lần) với xét nghiệm 6 tháng/lần.

Ba tỷ suất của mô hình được chọn thủ công để khớp số liệu công bố thời dolutegravir: tỷ suất bật lại ở một đoàn hệ miền Bắc (Matsumoto và c.s. 2025), tỷ lệ tái ức chế ở lần xét nghiệm lại (Sikombe và c.s. 2025) và mức kháng dolutegravir (World Health Organization 2024). Tỷ lệ ức chế virus toàn quốc không dùng để chọn tham số; Bài 3 dùng nó để đối chiếu. Đây không phải hiệu chuẩn hình thức (M08). Các tham số còn lại không có số liệu công khai phù hợp và được ghi rõ là giá trị minh họa trong bảng dưới.

flowchart LR
  S1["S1: Bậc 1 (TLD)<br/>ức chế virus"] -->|bật lại<br/>tùy tuân thủ| F1["F1: Bậc 1<br/>thất bại virus học"]
  F1 -->|tái ức chế<br/>nếu chưa kháng thuốc| S1
  F1 -.->|kháng thuốc mắc phải<br/>cờ của cá thể, cạnh tranh<br/>với tái ức chế và tử vong| F1
  F1 ==>|xét nghiệm phát hiện → tư vấn<br/>→ xét nghiệm lại vẫn cao| F2["F2: Bậc 2<br/>thất bại virus học"]
  F2 -->|ức chế| S2["S2: Bậc 2<br/>ức chế virus"]
  S2 -->|bật lại| F2
  S1 & F1 & S2 & F2 --> D["D: Tử vong"]

Mỗi cá thể mang thuộc tính cố định (tuổi khi vào mô hình, mức tuân thủ) và các biến “trí nhớ” thay đổi theo thời gian: cờ kháng thuốc, số tháng tích lũy sống với virus không ức chế, đồng hồ tư vấn tăng cường và lịch xét nghiệm kế tiếp. Chu kỳ một tháng, khung thời gian 70 năm (gần trọn đời), quan điểm hệ thống y tế (bên chi trả chính là BHYT). Chi phí tính bằng VNĐ năm 2025, lấy từ bảng giá dùng chung của khóa học: giá dịch vụ BHYT tại một bệnh viện hạng đặc biệt (Bộ Y tế 2024) và giá tham chiếu thuốc ARV của Global Fund quy đổi theo tỷ giá chính thức (The Global Fund 2026; Cục Thống kê (Bộ Tài chính) 2026). Tử vong nền lấy từ bảng sống Việt Nam của WHO (World Health Organization 2026). Kết cục là năm sống điều chỉnh theo chất lượng (QALY) và năm sống, chiết khấu 3%/năm; ngưỡng chi phí-hiệu quả cơ sở bằng 1 × GDP bình quân đầu người, 125,500,000 VNĐ/QALY (Cục Thống kê (Bộ Tài chính) 2026; Ong và c.s. 2026).

Bản đồ bài học

Bài Tiêu đề Nội dung chính Mã nguồn
1 Khi mô hình đoàn hệ không đủ Bất đẳng thức Jensen, cạn kiệt nhóm nguy cơ cao, bùng nổ số trạng thái, bất định ngẫu nhiên và bất định tham số code/l01-why-microsimulation.R
2 Mô phỏng vi mô bằng vòng lặp Thuộc tính cá thể, nguy cơ cạnh tranh với một số ngẫu nhiên, chuyển trạng thái phụ thuộc lịch sử, kiểm tra có thể thất bại, vết mô phỏng, dùng chung số ngẫu nhiên code/l02-microsim-loops.R
3 Vectơ hóa và đo hiệu năng Đảo thứ tự vòng lặp, thu về mô hình đoàn hệ, bench và bộ nhớ, Rprof, kiểm tra độ khớp và đối chiếu, phân tích cơ sở, giá xét nghiệm hòa vốn code/l03-vectorisation.R
4 Tính toán song song và Rcpp parLapply, mclapply, khối seed tách rời, C++ cho vòng lặp nóng, định luật Amdahl, máy không có trình biên dịch code/l04-parallel-rcpp.R, code/microsim.cpp
5 Sai số Monte Carlo, PSA và báo cáo Chọn số cá thể theo power, khi lợi ích sức khỏe dồn vào ít người, dùng chung số ngẫu nhiên, PSA lồng hai vòng, phân tách phương sai, CHEERS 2022 code/l05-mc-error-psa.R

Bảng tham số

Tham số riêng của module nằm trong params/hiv.csv. Cột “Nguồn” ghi trích dẫn đã kiểm chứng; “minh họa” nghĩa là giá trị minh họa, không có số liệu công khai phù hợp; “lựa chọn mô hình” là quyết định thiết kế (khung thời gian, khoảng tuổi). Các dòng tgt_* là mục tiêu kiểm tra và đối chiếu, chỉ dùng để so sánh đầu ra; dòng n_on_art chỉ dùng để quy kết quả lên cả nước. Với PSA (Bài 5), khoảng thấp-cao được hiểu là khoảng 95% và tham số phân phối được tính trong code bằng dist_params(), vì vậy cột dist_p1, dist_p2 để trống.

Tên Tham số Giá trị Thấp–cao Phân phối Đơn vị Nguồn
age_mean Tuổi trung bình của người đang điều trị ARV 38.2 fixed năm (Dat và c.s. 2018)
age_sd Độ lệch chuẩn của tuổi 7 fixed năm (Dat và c.s. 2022)
age_min Tuổi nhỏ nhất khi vào mô hình 18 fixed năm lựa chọn mô hình
age_max Tuổi lớn nhất khi vào mô hình 70 fixed năm lựa chọn mô hình
adh_mean Mức tuân thủ điều trị trung bình (tỷ lệ liều thuốc được uống) 0.9 0.85 – 0.95 beta tỷ lệ minh họa
adh_kappa Độ tập trung của phân phối Beta của tuân thủ 12 fixed không đơn vị minh họa
r_reb1 Tỷ suất bật lại tải lượng virus trên phác đồ bậc một (TLD) khi tuân thủ 100% 0.009 0.0045 – 0.018 gamma mỗi năm (Matsumoto và c.s. 2025)
hr_adh10 Tỷ số nguy cơ bật lại khi tuân thủ giảm 10 điểm phần trăm 2 1.5 – 2.7 lognormal tỷ số minh họa
r_resup1 Tỷ suất tái ức chế virus trên bậc một khi tuân thủ 100% (chưa kháng thuốc) 6.5 3.25 – 9.75 gamma mỗi năm (Sikombe và c.s. 2025)
hr_eac Hệ số tăng tái ức chế trong thời gian tư vấn tăng cường tuân thủ 2 1.2 – 3 lognormal tỷ số minh họa
eac_months Thời gian tư vấn tăng cường trước khi xét nghiệm lại 3 fixed tháng (Bộ Y tế 2021)
eac_sessions Số buổi tư vấn tăng cường tuân thủ trong một đợt 3 fixed buổi (Bộ Y tế 2021)
visits_per_year Số lượt khám và lĩnh thuốc ARV mỗi năm khi ổn định 4 fixed lượt/năm (Bộ Y tế 2021)
r_adr Tỷ suất xuất hiện kháng thuốc với phác đồ bậc một (kháng DTG) khi thất bại virus học 0.12 0.06 – 0.24 gamma mỗi năm (World Health Organization 2024)
r_reb2 Tỷ suất bật lại tải lượng virus trên phác đồ bậc hai khi tuân thủ 100% 0.03 0.015 – 0.06 gamma mỗi năm minh họa
r_resup2 Tỷ suất ức chế virus trên phác đồ bậc hai khi tuân thủ 100% 6 3 – 9 gamma mỗi năm minh họa
mu_supp Tỷ suất tử vong do HIV khi ức chế virus 0.004 0.002 – 0.008 gamma mỗi năm minh họa
mu_vir Tỷ suất tử vong do HIV khi thất bại virus học 0.04 0.02 – 0.08 gamma mỗi năm minh họa
vir_slope Mức tăng tương đối tử vong do HIV cho mỗi năm tích lũy sống với virus không ức chế 0.5 0.2 – 1 gamma mỗi năm tích lũy minh họa
horizon_years Khung thời gian mô phỏng 70 fixed năm lựa chọn mô hình
vl_soc Khoảng cách xét nghiệm tải lượng virus, chiến lược hiện hành 12 fixed tháng (Bộ Y tế 2021)
vl_alt Khoảng cách xét nghiệm tải lượng virus, chiến lược thay thế 6 fixed tháng (Pollack và c.s. 2019)
art2_usd_year Giá thuốc ARV bậc hai AZT + 3TC + LPV/r mỗi người-năm 279.6 181.44 – 377.76 gamma USD/người-năm (The Global Fund 2026)
vir_beddays_month Số ngày giường điều trị thêm mỗi tháng khi thất bại virus học (nhiễm trùng cơ hội, nhập viện) 1.5 0.75 – 3 gamma ngày giường/tháng minh họa
u_supp Giá trị thỏa dụng khi ức chế virus 0.8 0.78 – 0.82 beta thỏa dụng (Tran và c.s. 2018)
u_vir_dec Mức giảm thỏa dụng khi thất bại virus học 0.02 0.01 – 0.05 beta thỏa dụng minh họa
tgt_reb_rate Mục tiêu kiểm tra: tỷ suất bật lại trung bình lúc vào mô hình 0.022 0.022 – 0.032 fixed mỗi người-năm (Matsumoto và c.s. 2025)
tgt_resup_eac Mục tiêu kiểm tra: tỷ lệ tái ức chế ở lần xét nghiệm lại sau thất bại virus học, phác đồ bậc một (TLD) 0.85 0.839 – 0.861 fixed tỷ lệ (Sikombe và c.s. 2025)
tgt_res_f1 Mục tiêu kiểm tra: tỷ lệ kháng DTG trong nhóm không ức chế virus đang dùng phác đồ bậc một 0.0625 0.039 – 0.086 fixed tỷ lệ (World Health Organization 2024)
tgt_supp Đối chiếu: tỷ lệ người điều trị ARV có tải lượng dưới ngưỡng ức chế 0.96 0.96 – 0.984 fixed tỷ lệ (Báo Lào Cai 2025)
n_on_art Số người đang điều trị ARV tại Việt Nam (12/2024) 184,214 fixed người (Báo Lào Cai 2025)

Các giá trị dùng chung của khóa học không được chép vào params/hiv.csv; load_hiv_params() trong code/model.R đọc chúng từ kb/_shared mỗi lần chạy và giữ cố định trong PSA.

Tên Giá trị dùng chung Giá trị Đơn vị Nguồn
wtp Ngưỡng chi phí-hiệu quả cơ sở = 1 × GDP bình quân đầu người 125,500,000 VNĐ/QALY hoặc VNĐ/DALY tránh được (Ong và c.s. 2026)
wtp_low Ngưỡng thấp theo chi phí cơ hội = 0.5 × GDP bình quân đầu người 62,750,000 VNĐ/QALY hoặc VNĐ/DALY tránh được (Woods và c.s. 2016)
disc_rate Tỷ lệ chiết khấu hằng năm cho chi phí và hiệu quả 0.03 mỗi năm (Ong và c.s. 2026)
fx_vnd_per_usd Tỷ giá chính thức bình quân năm 2025 24,970 VNĐ/USD (Cục Thống kê (Bộ Tài chính) 2026)
art_tld_year Giá thuốc ARV TLD (tenofovir/lamivudin/dolutegravir) cho một người bệnh trong một năm 819,016 VNĐ/người bệnh-năm (The Global Fund 2026)
test_hiv_viral_load Giá xét nghiệm tải lượng HIV (hệ thống tự động) 979,700 VNĐ/xét nghiệm (Bộ Y tế 2024)
visit_outpatient Giá khám bệnh BHYT tại bệnh viện hạng đặc biệt (một lượt khám) 50,600 VNĐ/lượt khám (Bộ Y tế 2024)
bed_day_ward_infectious Giá ngày giường khoa Truyền nhiễm (nội khoa loại 1) tại bệnh viện hạng đặc biệt 327,100 VNĐ/ngày giường (Bộ Y tế 2024)

Tử vong nền: bảng sống Việt Nam năm 2023 của WHO, hai giới, qua vn_lifetable() (World Health Organization 2026).

Gói R và công cụ

  • R cơ bản cho cơ chế mô hình; ggplot2 cho hình; bench để đo hiệu năng (Hester và Vaughan 2025); parallel (có sẵn trong R) cho tính toán song song.
  • Rcpp cho bộ máy C++ (Eddelbuettel và François 2011). Cần trình biên dịch: trên macOS cài Xcode Command Line Tools (xcode-select --install; nếu Xcode đã cài nhưng chưa chấp nhận giấy phép, compile_microsim_cpp() tự dùng Command Line Tools cho lần biên dịch đó), trên Windows cài Rtools. Nếu chưa có, Bài 4 bỏ qua phần C++ kèm hướng dẫn cài đặt, còn Bài 5 tự chuyển sang bộ máy R (chậm hơn, cho cùng kết quả với cùng seed) và giảm số bộ tham số PSA.
  • Chạy mã của một bài: cd kb/m04-microsimulation && Rscript code/l03-vectorisation.R. Mỗi tệp chạy độc lập; hình được lưu vào figs/.

Tài liệu tham khảo