Bài 5: Xây cây quyết định từ đầu trong R

Từ số liệu công bố đến xác suất nhánh, giá trị kỳ vọng, ICER, phân tích độ nhạy một chiều và phân tích ngưỡng

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

Hai giờ sáng, ICU nhận một người bệnh sốc nhiễm khuẩn nghi do vi khuẩn Gram âm. Ở ICU Bệnh viện Trung ương Quân đội 108, cứ khoảng 3 người vào khoa thì có một người đã mang sẵn trong ruột vi khuẩn Gram âm kháng carbapenem (Toan và c.s. 2025). Bác sĩ trực phải chọn ngay: thêm colistin, mỗi ngày khoảng 3.64 triệu VNĐ (khoảng 15 ngày lương tối thiểu vùng I) và nguy cơ độc thận 36% (Bệnh viện Trung ương Thái Nguyên 2023; Eljaaly và c.s. 2021; Chính phủ 2025), hay chờ kháng sinh đồ vài ngày, hay gửi một xét nghiệm PCR giá 771,700 VNĐ (Bộ Y tế 2024). Mỗi lựa chọn mở ra một chuỗi khả năng: kháng hay không, điều trị sớm hay muộn, tổn thương thận hay không, sống hay chết. Bài này dạy bạn vẽ chuỗi đó thành một cây quyết định, tính nó trong R, rồi hỏi kết luận vững đến đâu.

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

Sau bài này, bạn có thể:

  1. Mã hóa một cây gồm nút cơ hội (chance node) và nút kết thúc (terminal node), với xác suất có điều kiện.
  2. Chuyển số liệu công bố thành xác suất nhánh và kiểm tra cây tái tạo đúng đầu vào.
  3. Tính giá trị kỳ vọng bằng tính ngược (rollback), rồi suy ra ICER và NMB.
  4. Làm OWSA, kịch bản và phân tích ngưỡng với câu hỏi đúng khi có nhiều chiến lược.
  5. Nhận ra khi nào cây không đủ.

Ý chính

  • Cây gồm nút quyết định (decision node), nút cơ hội và nút kết thúc. Xác suất ở nút cơ hội là xác suất có điều kiện; các nhánh của một nút cộng lại bằng 1.
  • Xác suất một đường đi là tích các xác suất nhánh; giá trị kỳ vọng là tổng (xác suất × giá trị nút kết thúc), và bằng kết quả tính ngược từ ngọn về gốc.
  • Chiến lược được mô tả bằng những gì nó đổi trên cây: hai xác suất và một chi phí ban đầu. Phủ CR (kháng carbapenem, carbapenem-resistant) cho mọi ca là “xét nghiệm” có độ nhạy 1, độ đặc hiệu 0. Kiểm tra trên đầu ra của cây: thường quy phải tái tạo tử vong và số ngày ICU quan sát.
  • Với từ ba chiến lược trở lên, OWSA và ngưỡng phải so với phương án tốt nhất còn lại, không chỉ với thường quy (Bài 3).
  • Khoảng thấp–cao có nghĩa khác nhau (khoảng tin cậy, phân tán giữa người bệnh, loại thuốc hay loại xét nghiệm khác); nghi ngờ cấu trúc đi vào kịch bản và ngưỡng. Trong ví dụ, không tham số đơn lẻ nào làm đổi kết luận; chỉ một kịch bản gộp nhiều giả định bất lợi mới làm đổi.
  • Cây hợp với vấn đề ngắn hạn, không lặp lại; sau đó cần mô hình Markov.

Nội dung

Cấu trúc cây

Mỗi chiến lược mở ra cùng một cây con gồm bốn tầng nút cơ hội, mỗi nút kết thúc mang số ngày ICU, chi phí, năm sống, QALY và DALY (Detsky và c.s. 1997; Roberts và c.s. 2012):

flowchart LR
  D["Chiến lược<br/>(nút quyết định)"] --> C1(("Kháng<br/>carbapenem?"))
  C1 -->|"có: p_cr"| C2(("Điều trị phù<br/>hợp sớm?"))
  C1 -->|"không: 1 − p_cr"| C3(("Bị phủ CR<br/>không cần thiết?"))
  C2 -->|"sớm"| M1(("Tử vong?"))
  C2 -->|"muộn"| M2(("Tử vong?"))
  C3 -->|"có"| A(("Tổn thương<br/>thận cấp?"))
  C3 -->|"không"| M5(("Tử vong?"))
  A -->|"có"| M3(("Tử vong?"))
  A -->|"không"| M4(("Tử vong?"))
  M1 --> E1(["Sống / Chết"])
  M2 --> E2(["Sống / Chết"])
  M3 --> E3(["Sống / Chết"])
  M4 --> E4(["Sống / Chết"])
  M5 --> E5(["Sống / Chết"])

flowchart LR
  D["Chiến lược<br/>(nút quyết định)"] --> C1(("Kháng<br/>carbapenem?"))
  C1 -->|"có: p_cr"| C2(("Điều trị phù<br/>hợp sớm?"))
  C1 -->|"không: 1 − p_cr"| C3(("Bị phủ CR<br/>không cần thiết?"))
  C2 -->|"sớm"| M1(("Tử vong?"))
  C2 -->|"muộn"| M2(("Tử vong?"))
  C3 -->|"có"| A(("Tổn thương<br/>thận cấp?"))
  C3 -->|"không"| M5(("Tử vong?"))
  A -->|"có"| M3(("Tử vong?"))
  A -->|"không"| M4(("Tử vong?"))
  M1 --> E1(["Sống / Chết"])
  M2 --> E2(["Sống / Chết"])
  M3 --> E3(["Sống / Chết"])
  M4 --> E4(["Sống / Chết"])
  M5 --> E5(["Sống / Chết"])

Từ số liệu công bố đến xác suất nhánh

prepare_inputs() thực hiện bốn phép chuyển đổi:

  1. Tách tử vong theo tình trạng kháng. Tử vong thô 31.7% (Enterobacterales kháng) so với 20.2% (nhạy), tỷ số 1.57 (Larsson và c.s. 2025), áp cho mọi chủng kháng. Chia tử vong ICU 40.1% (Do và c.s. 2021) thành hai nhóm giữ đúng tỷ số và trung bình: 34.2% (không kháng), 53.7% (kháng, điều trị muộn).
  2. Áp tỷ số chênh 1.60 của điều trị không phù hợp (Paul và c.s. 2010): chia odds (số chênh, \(p/(1-p)\)), được tử vong 42.0% khi điều trị sớm. Hướng tác động hợp lý (Kumar và c.s. 2006), nhưng ta áp OR này cho cả khoảng chờ kháng sinh đồ.
  3. Độc thận dư thừa từ colistin: 36.2%, nguy cơ tương đối 2.40 (Eljaaly và c.s. 2021), nhân hệ số 0.5 cho đợt ngắn (giá trị minh họa), còn 10.6%. Tổn thương thận cấp nhân odds tử vong với 2.45, tỷ số chênh gộp ở người bệnh viêm phổi nặng (Xie và c.s. 2024); đây là liên quan thống kê nên có thể phóng đại tác hại thật của vài ngày colistin. Tử vong của người không kháng bị độc thận vì vậy là 56.0% thay vì 34.2%.
  4. Số ngày ICU. Nguồn báo cáo trung bình 14.3 ngày cho mọi người bệnh, cả sống và tử vong (Hung và c.s. 2022). Ta suy ra 16.7 ngày cho người sống sót để thường quy tái tạo đúng trung bình.

Tỷ lệ kháng 30.2% tính theo chủng Gram âm phân lập ở 324 bệnh viện tự nguyện, không riêng ICU (Larsson và c.s. 2025), và khác nhau nhiều giữa các loài (Vu và c.s. 2021). Xét nghiệm đo gì? Độ nhạy 95% và độ đặc hiệu 99% là của nhóm nghiên cứu trên phết trực tràng, phát hiện gen carbapenemase ở vi khuẩn định cư, không phải tác nhân gây nhiễm khuẩn huyết (Bai và c.s. 2021). Nhóm chủng phân lập có độ nhạy cao hơn nhưng cần cấy dương tính trước, nên chỉ rút ngắn khoảng từ cấy đến kháng sinh đồ. Dùng các số này cho nhiễm khuẩn huyết do chủng kháng là giả định lạc quan; cov_cr (giá trị minh họa, cơ sở bằng 1) nhân với độ nhạy để thử hệ quả.

Mỗi chiến lược thay đổi gì trên cây

Chiến lược P(điều trị sớm | kháng) P(phủ CR không cần | không kháng)
Thường quy 0 0
Xét nghiệm nhanh cov_cr × Se × xác suất hành động (1 − Sp) × xác suất hành động
Xét nghiệm nhanh + AMS như trên, xác suất hành động cao hơn như trên
Phủ CR cho mọi ca 1 1

Xác suất bác sĩ hành động là 95% có AMS, 60% không có (giá trị minh họa).

Kết quả, OWSA, ngưỡng và kịch bản

Chi phí kỳ vọng của xét nghiệm nhanh kèm AMS, tính ngược hay theo đường đi, là 100.232 triệu VNĐ. ICER so với thường quy là 1.76 triệu VNĐ/QALY, NMB gia tăng 48.53 triệu VNĐ ở ngưỡng của khóa học. Với bốn chiến lược, câu hỏi đúng là so với phương án tốt nhất còn lại: ở cơ sở đó là Xét nghiệm nhanh (bị lấn át mở rộng nhưng NMB vẫn cao hơn thường quy), và biên còn 18.03 triệu VNĐ.

OWSA đặt từng tham số ở giá trị thấp rồi cao (Briggs và c.s. 2012) và đo biên đó. Cột kind cho biết khoảng thấp–cao là gì: khoảng tin cậy 95%, khoảng tứ phân vị (phân tán giữa người bệnh, không phải bất định của trung bình; xem M03 Bài 1), hai cách định giá, loại xét nghiệm hay loại thuốc khác, hay giá trị minh họa. Hai tham số đứng đầu là act_no_ams và cov_cr (Hình 12.1), đều là giá trị minh họa. act_no_ams đứng đầu vì nó quyết định đối thủ mạnh nhất, xét nghiệm nhanh không có AMS, tốt đến đâu. Không tham số nào làm đổi phương án tối ưu trong khoảng của nó: biên nhỏ nhất là 1.72 triệu VNĐ, khi cov_cr ở giá trị thấp. Giá xét nghiệm c_test không có trong biểu đồ vì biên không đổi trong khoảng thấp–cao: đối thủ mạnh nhất, xét nghiệm nhanh không có AMS, trả đúng khoản đó. Giá chỉ quan trọng khi đối thủ mạnh nhất đổi, như phần phân tích ngưỡng ngay sau đây.

Phân tích ngưỡng (threshold analysis) tìm giá trị tham số làm biên bằng 0. Biên cắt 0 hai lần theo tỷ lệ kháng (Hình 12.2): xét nghiệm nhanh kèm AMS tối ưu khi tỷ lệ kháng từ 0.9% đến 71.3%. Tỷ lệ kháng của cả E. coli (11%), K. pneumoniae (37%) và A. baumannii (64%) đều nằm trong khoảng (Larsson và c.s. 2025). Giá xét nghiệm tối đa là 49.30 triệu VNĐ nếu chỉ so với thường quy, nhưng 27.41 triệu khi so với cả bốn chiến lược, vì khi giá tăng Phủ CR cho mọi ca vượt lên. Cả hai đều cao hơn nhiều so với giá BHYT 771,700 VNĐ.

Kịch bản (chunk l05-scenarios) dành cho nghi ngờ mà khoảng thấp–cao không diễn đạt được. Khi A. baumannii chiếm ưu thế và xét nghiệm chỉ bắt một nửa ca (ví dụ khi cơ chế kháng nằm ngoài bộ gen mà xét nghiệm phát hiện), tối ưu là Phủ CR cho mọi ca. Khi xét nghiệm bắt một nửa ca và hiệu quả điều trị sớm chỉ bằng nửa (giả định minh họa), ICER so với thường quy lên 5.68 triệu VNĐ/QALY và biên giảm còn 4.44 triệu VNĐ, nhưng xét nghiệm nhanh kèm AMS vẫn tối ưu. Kết luận của cây vì vậy vững với từng tham số riêng lẻ, nhưng vẫn có điều kiện về loài vi khuẩn và về những gì xét nghiệm thực sự phát hiện.

Khi nào cây quyết định không đủ

Cây giả định sự kiện xảy ra ngắn hạn và không lặp lại; giả định gãy khi có tái nhiễm, nguy cơ đổi theo thời gian, hoặc khung thời gian dài. Một sự kiện có/không mỗi tháng cần \(2^n\) đường đi sau \(n\) tháng: 4,096 sau một năm. Mô hình Markov thay các nhánh lặp bằng một tập trạng thái nhỏ và một ma trận chuyển dùng lại mỗi chu kỳ (M02) (Sonnenberg và Beck 1993; Barton, Bryan, và Robinson 2004). Khi người bệnh lây cho nhau, cần mô hình lây truyền động (M05).

Thực hành

Nạp hàm và tham số.

if (!isTRUE(l10n_info()[["UTF-8"]])) invisible(Sys.setlocale("LC_CTYPE", "en_US.UTF-8"))
source("../_shared/R/helpers.R")
source("../_shared/R/vn_data.R")   # course-wide values: threshold, life table, CPI, prices
source("code/model.R")
library(ggplot2)
p <- load_params()                   # module parameters + shared values, costs in 2025 VND
inp <- prepare_inputs(p)

Xác suất suy ra, kèm kiểm tra trên đầu ra của cây.

# From published summaries to branch probabilities (see prepare_inputs() in model.R)
derived <- c(rr_cr = inp$rr_cr, m_noncr = inp$m_noncr, m_cr_late = inp$m_cr_late,
             m_cr_early = inp$m_cr_early, p_aki_excess = inp$p_aki_excess,
             los_survivor_days = inp$los_surv)
round(derived, 4)
            rr_cr           m_noncr         m_cr_late        m_cr_early 
           1.5693            0.3422            0.5370            0.4202 
     p_aki_excess los_survivor_days 
           0.1056           16.6722 
# Checks on the TREE OUTPUT, not on the derivation formulas: usual care must reproduce
# the observed mortality and mean ICU days, the CR branch must carry mass p_cr, and the
# odds ratio between late and early therapy must be the published one
paths_usual <- tree_paths(p, "usual")
tree_check <- c(mortality_tree = sum(paths_usual$prob[paths_usual$dead]),
                mortality_data = p$m_sepsis_icu,
                icu_days_tree = sum(paths_usual$prob * paths_usual$los),
                icu_days_data = p$los_icu,
                cr_mass_tree = sum(paths_usual$prob[paths_usual$cr == "CR"]),
                p_cr = p$p_cr,
                or_from_tree = odds(inp$m_cr_late) / odds(inp$m_cr_early),
                or_input = p$or_inapprop)
round(tree_check, 4)
mortality_tree mortality_data  icu_days_tree  icu_days_data   cr_mass_tree 
         0.401          0.401         14.300         14.300          0.302 
          p_cr   or_from_tree       or_input 
         0.302          1.600          1.600 
stopifnot(abs(tree_check[["mortality_tree"]] - tree_check[["mortality_data"]]) < 1e-12,
          abs(tree_check[["icu_days_tree"]] - tree_check[["icu_days_data"]]) < 1e-9,
          abs(tree_check[["cr_mass_tree"]] - tree_check[["p_cr"]]) < 1e-12,
          abs(tree_check[["or_from_tree"]] - tree_check[["or_input"]]) < 1e-12)

Bảng chiến lược: mỗi hàng chỉ có hai xác suất và một chi phí ban đầu.

# Each strategy changes only two branch probabilities and the up-front cost
st <- strategy_table(p)
st_print <- st
st_print[, c("p_early_cr", "p_exposed_noncr")] <- round(st_print[, c("p_early_cr", "p_exposed_noncr")], 4)
st_print
   strategy                  label p_early_cr p_exposed_noncr cost_upfront
1     usual             Thường quy     0.0000          0.0000          0.0
2       rdt       Xét nghiệm nhanh     0.5700          0.0060     771700.0
3   rdt_ams Xét nghiệm nhanh + AMS     0.9025          0.0095     985801.9
4 cover_all      Phủ CR cho mọi ca     1.0000          1.0000          0.0

Mười đường đi của xét nghiệm nhanh kèm AMS; prob cộng bằng 1.

# Every root-to-leaf path for "rapid test + AMS": conditional probabilities at the
# four chance nodes, the path probability (their product) and the payoffs
paths_ams <- tree_paths(p, "rdt_ams")
show_cols <- c("cr", "therapy", "aki", "dead", "prob", "los", "cost", "qaly")
paths_print <- paths_ams[, show_cols]
paths_print$prob <- round(paths_print$prob, 5)
paths_print$los <- round(paths_print$los, 1)
paths_print$cost <- round(paths_print$cost / 1e6, 2)
paths_print$qaly <- round(paths_print$qaly, 2)
paths_print
      cr  therapy   aki  dead    prob  los   cost  qaly
1     CR    early FALSE  TRUE 0.11453 10.0  99.34  0.00
2     CR    early FALSE FALSE 0.15802 16.7 140.64 12.39
3     CR     late FALSE  TRUE 0.01581 11.0 105.53  0.00
4     CR     late FALSE FALSE 0.01363 17.7 146.83 12.39
5  nonCR  exposed  TRUE  TRUE 0.00039 13.0  92.43  0.00
6  nonCR  exposed  TRUE FALSE 0.00031 19.7 133.72 12.39
7  nonCR  exposed FALSE  TRUE 0.00203 10.0  73.85  0.00
8  nonCR  exposed FALSE FALSE 0.00390 16.7 115.14 12.39
9  nonCR standard FALSE  TRUE 0.23657 10.0  62.92  0.00
10 nonCR standard FALSE FALSE 0.45480 16.7 104.22 12.39
# Structural checks for every strategy: valid probabilities, leaves sum to 1, costs >= 0
for (s in STRATEGIES) check_tree(tree_paths(p, s))

Tính ngược từng nút; kết quả phải trùng tổng theo đường đi.

# Rollback ("folding back"): replace each chance node by its expected value,
# from the leaves to the root. It must equal the sum over paths of prob x payoff.
pa <- paths_ams
ev_leaf <- function(rows) sum(pa$p_node4[rows] * pa$cost[rows])   # death node
ev_early <- ev_leaf(1:2)
ev_late <- ev_leaf(3:4)
ev_aki <- ev_leaf(5:6)
ev_noaki <- ev_leaf(7:8)
ev_notexp <- ev_leaf(9:10)
ev_exposed <- inp$p_aki_excess * ev_aki + (1 - inp$p_aki_excess) * ev_noaki          # AKI node
s_ams <- st[st$strategy == "rdt_ams", ]
ev_cr <- s_ams$p_early_cr * ev_early + (1 - s_ams$p_early_cr) * ev_late             # therapy node, CR
ev_noncr <- s_ams$p_exposed_noncr * ev_exposed + (1 - s_ams$p_exposed_noncr) * ev_notexp
ev_root <- p$p_cr * ev_cr + (1 - p$p_cr) * ev_noncr                                 # CR node
c(rollback = ev_root, path_sum = sum(pa$prob * pa$cost)) / 1e6
rollback path_sum 
100.2321 100.2321 
stopifnot(abs(ev_root - sum(pa$prob * pa$cost)) < 1e-6)

Kết quả kỳ vọng, ICER, NMB gia tăng và biên so với đối thủ mạnh nhất.

# Expected results for all strategies, then the frontier
res <- run_tree(p)
res_print <- data.frame(strategy = res$strategy,
                        cost_million = res$cost / 1e6, deaths_per_1000 = 1000 * res$death,
                        qaly = res$qaly, daly = res$daly, nmb_million = res$nmb / 1e6)
res_print[, -1] <- round(res_print[, -1], 3)
res_print
   strategy cost_million deaths_per_1000  qaly  daly nmb_million
1     usual       99.541         401.000 7.419 7.698     831.546
2       rdt      100.127         381.000 7.667 7.314     862.048
3   rdt_ams      100.232         369.333 7.811 7.090     880.078
4 cover_all      107.458         381.820 7.657 7.330     853.443
fr <- efficiency_frontier(res$strategy, res$cost, res$qaly)
icer_main <- fr$icer[fr$strategy == "rdt_ams"]
inmb_base <- inmb(p)                 # incremental NMB versus usual care only
competitor_base <- best_other(p)     # strongest rival of rapid test + AMS when ALL strategies compete
margin_base <- nmb_margin(p)         # NMB advantage over that rival
round(c(icer_million = icer_main / 1e6, inmb_vs_usual_million = inmb_base / 1e6,
        margin_vs_best_other_million = margin_base / 1e6), 3)
                icer_million        inmb_vs_usual_million 
                       1.762                       48.532 
margin_vs_best_other_million 
                      18.029 
competitor_base
[1] "rdt"
# The prose says the closest rival is off the frontier (dominated or extendedly dominated),
# yet has a higher NMB than usual care
stopifnot(competitor_base != "usual", fr$status[fr$strategy == competitor_base] != "frontier")

OWSA trên biên, kèm chiến lược tối ưu ở mỗi đầu; dấu của biên phải khớp với chiến lược tối ưu.

# One-way sensitivity analysis. With 4 strategies the outcome must compare rapid test + AMS
# with the BEST of the other three (nmb_margin), not with usual care alone. Besides the
# margin we record which strategy is optimal at each bound.
tab <- attr(p, "table")
# Every input with a low and a high value, except the threshold itself (Lesson 4 covers it)
owsa_names <- setdiff(tab$name[!is.na(tab$low) & !is.na(tab$high)], "wtp")
ow <- owsa(p, owsa_names, outcome = nmb_margin)
# What each low-high range means (it decides how to read the tornado)
special_kind <- c(c_icu_nonbed_day = "IQR", bhyt_share = "IQR", c_nonmed = "định giá", disc_rate = "quy ước",
                  c_test = "loại xét nghiệm", c_cr_drug_day = "loại thuốc", c_ams = "giả định", u_post = "nghiên cứu khác")
illustrative <- as.logical(tab$illustrative[match(ow$name, tab$name)])
ow$kind <- ifelse(ow$name %in% names(special_kind), special_kind[ow$name],
                  ifelse(illustrative, "minh họa", "KTC 95%"))
# Parameters that cannot change this comparison have a range of zero
ow$name[ow$range < 1e-6]
[1] "dur_cr_course" "bhyt_share"    "c_nonmed"      "dw_sepsis"    
[5] "c_test"       
# The prose explains why the test price is absent from the tornado: its rival (rapid test without
# AMS) pays the same price, so the margin does not move
stopifnot("c_test" %in% ow$name[ow$range < 1e-6], competitor_base == "rdt")
ow_top <- ow[ow$range >= 1e-6, ][1:10, ]
ow_print <- ow_top
ow_print[, c("out_low", "out_high", "range")] <- round(ow_print[, c("out_low", "out_high", "range")] / 1e6, 2)
ow_print$low <- format(ow_print$low, big.mark = ",", scientific = FALSE, trim = TRUE, drop0trailing = TRUE)
ow_print$high <- format(ow_print$high, big.mark = ",", scientific = FALSE, trim = TRUE, drop0trailing = TRUE)
ow_print[, c("name", "kind", "low", "high", "out_low", "out_high", "best_low", "best_high")]
               name            kind  low high out_low out_high best_low
7        act_no_ams        minh họa  0.4  0.8   26.63     7.60  rdt_ams
6            cov_cr        minh họa  0.5    1    1.72    18.03  rdt_ams
12     or_aki_death         KTC 95% 1.14 5.26    5.86    17.96  rdt_ams
3       or_inapprop         KTC 95% 1.37 1.86   12.22    23.53  rdt_ams
26        disc_rate         quy ước    0 0.05   24.77    14.95  rdt_ams
24           u_post nghiên cứu khác 0.56 0.89   13.44    21.41  rdt_ams
8           act_ams        minh họa 0.85    1   12.82    20.64  rdt_ams
23   le_post_factor        minh họa  0.6    1   14.56    20.96  rdt_ams
11 aki_short_factor        minh họa 0.25    1   13.91    17.94  rdt_ams
10  rr_aki_colistin         KTC 95% 1.47 3.91   15.13    18.00  rdt_ams
   best_high
7    rdt_ams
6    rdt_ams
12   rdt_ams
3    rdt_ams
26   rdt_ams
24   rdt_ams
8    rdt_ams
23   rdt_ams
11   rdt_ams
10   rdt_ams
# Which single parameters change the optimal strategy at one of their bounds?
flip_names <- ow$name[ow$best_low != "rdt_ams" | ow$best_high != "rdt_ams"]
flip_names
character(0)
# No single parameter changes the optimal strategy within its range (the prose says so);
# the smallest margin reached at any bound is still positive
margin_min <- min(c(ow$out_low, ow$out_high))
ow_min <- ow$name[which.min(pmin(ow$out_low, ow$out_high))]
c(smallest_margin_million = round(margin_min / 1e6, 2), at = ow_min)
smallest_margin_million                      at 
                 "1.72"                "cov_cr" 
stopifnot(length(flip_names) == 0, margin_min > 0, ow$out_low[ow$name == ow_min] == margin_min)
# The prose explains the top of the tornado: act_no_ams first, and the top two are illustrative
stopifnot(ow_top$name[1] == "act_no_ams", all(ow_top$kind[1:2] == "minh họa"))
# The sign of the margin and the identity of the optimal strategy must agree
stopifnot(all((ow$out_low > 0) == (ow$best_low == "rdt_ams")),
          all((ow$out_high > 0) == (ow$best_high == "rdt_ams")))

ow_top$label <- paste0(ow_top$name, " (", ow_top$kind, ")")
ow_top$label <- factor(ow_top$label, levels = rev(ow_top$label))
fig_tornado <- ggplot(ow_top) +
  geom_segment(aes(x = label, xend = label, y = margin_base / 1e6, yend = out_low / 1e6,
                   colour = "Giá trị thấp"), linewidth = 6) +
  geom_segment(aes(x = label, xend = label, y = margin_base / 1e6, yend = out_high / 1e6,
                   colour = "Giá trị cao"), linewidth = 6) +
  geom_hline(yintercept = margin_base / 1e6, colour = "grey30") +
  geom_hline(yintercept = 0, linetype = "22", colour = pal_course[2]) +
  coord_flip() +
  scale_colour_manual(values = c("Giá trị thấp" = pal_course[5], "Giá trị cao" = pal_course[1]), name = NULL) +
  labs(x = NULL, y = "Biên NMB so với phương án tốt nhất còn lại (triệu VNĐ)",
       title = "Biểu đồ tornado: tham số nào quyết định kết quả?",
       subtitle = "Dưới đường đứt (< 0): Xét nghiệm nhanh + AMS không còn là phương án tối ưu") +
  theme_course()
save_fig(fig_tornado, "l05-tornado.png", height = 5)
Hình 1: Tornado: biên NMB của xét nghiệm nhanh kèm AMS khi từng tham số ở giá trị thấp hoặc cao. Dưới đường đứt, nó không còn tối ưu.

Ngưỡng cho tỷ lệ kháng và giá xét nghiệm; ngay ngoài ngưỡng, chiến lược tối ưu phải đã đổi.

# Threshold analysis on the same outcome. The margin crosses zero twice along p_cr:
# at low resistance usual care wins, at high resistance cover-all wins.
pcr_low <- threshold_value(p, "p_cr", lower = 0.001, upper = p$p_cr)
pcr_high <- threshold_value(p, "p_cr", lower = p$p_cr, upper = 0.99)
c(p_cr_low = pcr_low, p_cr_high = pcr_high)
   p_cr_low   p_cr_high 
0.008869681 0.713361613 
# At each crossing the margin is zero, and just outside the optimal strategy has changed
stopifnot(abs(nmb_margin(modifyList(p, list(p_cr = pcr_low)))) < 1,
          abs(nmb_margin(modifyList(p, list(p_cr = pcr_high)))) < 1,
          best_at(modifyList(p, list(p_cr = pcr_low * 0.9))) == "usual",
          best_at(modifyList(p, list(p_cr = pcr_high * 1.1))) == "cover_all")

# All three species prevalences lie between the two thresholds (the prose says so)
stopifnot(pcr_low < p$p_cr_ecoli, p$p_cr_abau < pcr_high)

# Maximum test price: against usual care only, and against the best of all strategies
price_pair <- threshold_value(p, "c_test", lower = 0, upper = 100 * p$c_test, outcome = inmb)
price_full <- threshold_value(p, "c_test", lower = 0, upper = 100 * p$c_test)
competitor_above <- best_other(modifyList(p, list(c_test = 1.01 * price_full)))
round(c(base_price = p$c_test, max_price_vs_usual_only = price_pair,
        max_price_vs_all = price_full) / 1e6, 3)
             base_price max_price_vs_usual_only        max_price_vs_all 
                  0.772                  49.304                  27.407 
competitor_above
[1] "cover_all"
# Against all strategies the maximum price is lower, yet still far above the BHYT price
stopifnot(price_full < price_pair, price_full > 5 * p$c_test, competitor_above != "usual")

# Incremental NMB of each alternative versus usual care along the prevalence of resistance
pcr_grid <- seq(0.01, 0.80, by = 0.01)   # wide enough to show both crossings
inmb_grid <- NULL
for (x in pcr_grid) {
  r_x <- run_tree(modifyList(p, list(p_cr = x)))
  inmb_grid <- rbind(inmb_grid, data.frame(p_cr = x, strategy = r_x$label[-1],
                                           inmb = r_x$nmb[-1] - r_x$nmb[1]))
}
inmb_grid$strategy <- factor(inmb_grid$strategy, levels = STRATEGY_LABELS[-1])
fig_threshold <- ggplot(inmb_grid, aes(100 * p_cr, inmb / 1e6, colour = strategy)) +
  geom_hline(yintercept = 0, colour = "grey40") +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = 100 * c(pcr_low, pcr_high), linetype = "22", colour = "grey40") +
  geom_vline(xintercept = 100 * p$p_cr, linetype = "dotted", colour = "grey40") +
  scale_colour_manual(values = pal_course[c(4, 1, 2)], name = NULL) +
  labs(x = "Tỷ lệ kháng carbapenem trong chủng Gram âm (%)",
       y = "NMB gia tăng so với thường quy (triệu VNĐ)",
       title = "Phương án tối ưu theo tỷ lệ kháng carbapenem",
       subtitle = "Đường cao nhất là tối ưu; đường đứt: hai ngưỡng đổi quyết định; đường chấm: phân tích cơ sở") +
  theme_course()
save_fig(fig_threshold, "l05-threshold.png")
Hình 2: NMB gia tăng so với thường quy theo tỷ lệ kháng. Đường cao nhất là tối ưu; hai đường đứt là hai ngưỡng.

Kịch bản theo loài và giả định về xét nghiệm; hàng cơ sở khớp kết quả cơ sở.

# Scenario analysis: structural doubts that a low-high range cannot express. Settings
# marked (*) are illustrative scenario choices, not table parameters.
cov_low <- param_bounds(p, "cov_cr")[["low"]]
scenarios <- list(
  list(name = "Phân tích cơ sở", change = list()),
  list(name = "E. coli chiếm ưu thế", change = list(p_cr = p$p_cr_ecoli)),
  list(name = "K. pneumoniae chiếm ưu thế", change = list(p_cr = p$p_cr_kpn)),
  list(name = "A. baumannii chiếm ưu thế, xét nghiệm bắt được một nửa ca (*)", change = list(p_cr = p$p_cr_abau, cov_cr = cov_low)),
  list(name = "Hiệu quả điều trị sớm bằng nửa (OR^0.5) (*)", change = list(or_inapprop = sqrt(p$or_inapprop))),
  list(name = "Bi quan: xét nghiệm bắt một nửa ca và hiệu quả bằng nửa (*)", change = list(cov_cr = cov_low, or_inapprop = sqrt(p$or_inapprop)))
)
scen_tab <- data.frame(scenario = sapply(scenarios, function(x) x$name), p_cr = NA_real_,
                       best = NA_character_, frontier = NA_character_,
                       icer_vs_usual_million = NA_real_, margin_million = NA_real_)
for (k in seq_along(scenarios)) {
  q <- modifyList(p, scenarios[[k]]$change)
  r_q <- run_tree(q)
  fr_q <- efficiency_frontier(r_q$strategy, r_q$cost, r_q$qaly)
  scen_tab$p_cr[k] <- q$p_cr
  scen_tab$best[k] <- STRATEGY_LABELS[[best_at(q)]]
  scen_tab$frontier[k] <- paste(STRATEGY_LABELS[fr_q$strategy[fr_q$status == "frontier"]], collapse = " > ")
  scen_tab$icer_vs_usual_million[k] <- (r_q["rdt_ams", "cost"] - r_q["usual", "cost"]) /
    (r_q["rdt_ams", "qaly"] - r_q["usual", "qaly"]) / 1e6
  scen_tab$margin_million[k] <- nmb_margin(q) / 1e6
}
# The base-case row must reproduce the base results (checked before rounding for display)
stopifnot(abs(scen_tab$margin_million[1] - margin_base / 1e6) < 1e-9)
# Prose: the A. baumannii scenario makes cover-all optimal; in the pessimistic scenario rapid
# test + AMS stays optimal with a smaller margin
n_frontier_pess <- lengths(strsplit(scen_tab$frontier[6], " > "))
stopifnot(scen_tab$best[4] == STRATEGY_LABELS[["cover_all"]],
          scen_tab$best[6] == STRATEGY_LABELS[["rdt_ams"]],
          scen_tab$margin_million[6] < scen_tab$margin_million[1],
          sum(scen_tab$best != STRATEGY_LABELS[["rdt_ams"]]) == 1,   # only one scenario flips
          scen_tab$icer_vs_usual_million[6] > scen_tab$icer_vs_usual_million[1])   # pessimistic ICER is higher
scen_print <- scen_tab
scen_print[, c("p_cr", "icer_vs_usual_million", "margin_million")] <-
  round(scen_print[, c("p_cr", "icer_vs_usual_million", "margin_million")], 2)
scen_print
                                                       scenario p_cr
1                                               Phân tích cơ sở 0.30
2                                          E. coli chiếm ưu thế 0.11
3                                    K. pneumoniae chiếm ưu thế 0.37
4 A. baumannii chiếm ưu thế, xét nghiệm bắt được một nửa ca (*) 0.64
5                   Hiệu quả điều trị sớm bằng nửa (OR^0.5) (*) 0.30
6   Bi quan: xét nghiệm bắt một nửa ca và hiệu quả bằng nửa (*) 0.30
                    best
1 Xét nghiệm nhanh + AMS
2 Xét nghiệm nhanh + AMS
3 Xét nghiệm nhanh + AMS
4      Phủ CR cho mọi ca
5 Xét nghiệm nhanh + AMS
6 Xét nghiệm nhanh + AMS
                                                 frontier icer_vs_usual_million
1                     Thường quy > Xét nghiệm nhanh + AMS                  1.76
2                     Thường quy > Xét nghiệm nhanh + AMS                  6.79
3                     Thường quy > Xét nghiệm nhanh + AMS                  1.22
4 Thường quy > Xét nghiệm nhanh + AMS > Phủ CR cho mọi ca                  1.37
5                     Thường quy > Xét nghiệm nhanh + AMS                  0.19
6                     Thường quy > Xét nghiệm nhanh + AMS                  5.68
  margin_million
1          18.03
2           6.30
3          22.09
4         -46.45
5           9.21
6           4.44

Số đường đi khi có sự kiện lặp lại hằng tháng, so với số trạng thái Markov.

# When trees break: one yes/no event per month gives 2^n paths after n months,
# while a Markov model keeps the same few states every cycle (Module 2)
months <- c(1, 6, 12, 24, 60)
explosion <- data.frame(months = months, tree_paths = fmt_num(2^months), markov_states = 3)
explosion
  months                tree_paths markov_states
1      1                         2             3
2      6                        64             3
3     12                     4,096             3
4     24                16,777,216             3
5     60 1,152,921,504,606,846,976             3
Mẹo thực hành
  • Vẽ cây trên giấy trước khi viết code: mỗi tầng nút cơ hội là một cột xác suất.
  • Nạp mọi giá trị dùng chung (ngưỡng, giá BHYT, giá thuốc, bảng sống) qua một hàm duy nhất như load_params(); khi bảng giá mới ra, chỉ cần sửa một tệp chung và chạy lại.
  • Ghi năm giá cạnh mỗi chi phí (giá thầu 2023, giá BHYT 2025) và quy về cùng một năm trước khi cộng.
Lỗi thường gặp
  • Xác suất không có điều kiện, ví dụ dùng tử vong chung cho nhánh kháng và không kháng.
  • OWSA và ngưỡng chỉ so với thường quy khi có từ ba chiến lược: kết luận “robust” có thể sai.
  • Chia xác suất cho tỷ số chênh (Bài tập 3).
  • Quên người được cứu sống tiêu tốn thêm nguồn lực: người sống nằm ICU lâu hơn.

Bài tập

Bài 1 (dễ). Viết hai kiểm thử hộp đen. (i) Xét nghiệm và AMS miễn phí, không bác sĩ nào hành động: kết quả phải bằng thường quy. (ii) Nếu tỷ số chênh bằng 1, chênh lệch tử vong chỉ đến từ độc thận do colistin: tính bằng tay và so với cây.

Ý (i): mọi chênh lệch bằng 0. Ý (ii): cây cho 152.7 ca tử vong thêm trên một triệu người bệnh, bằng phép tính tay. Loại kiểm tra này bắt được nhiều lỗi chỉ số.

# Exercise 1: black-box tests. (i) A free test that nobody acts on equals usual care.
q <- modifyList(p, list(act_ams = 0, c_test = 0, c_ams = 0))
r_q <- run_tree(q)
d_q <- incremental(r_q, "rdt_ams", "usual")
round(d_q[c("cost", "qaly", "death")], 10)
 cost  qaly death 
    0     0     0 
stopifnot(all(abs(d_q[c("cost", "qaly", "death")]) < 1e-8))
# (ii) With odds ratio 1, early therapy brings no benefit, so any difference in deaths
# can only come from the excess kidney injury of unnecessary colistin; compute it by hand
q1 <- modifyList(p, list(or_inapprop = 1))
d1 <- incremental(run_tree(q1), "rdt_ams", "usual")
by_hand <- (1 - p$p_cr) * strategy_table(q1)$p_exposed_noncr[3] * inp$p_aki_excess * (inp$m_noncr_aki - inp$m_noncr)
round(1e6 * c(tree = d1[["death"]], by_hand = by_hand), 2)   # extra deaths per million patients
   tree by_hand 
 152.73  152.73 
stopifnot(abs(d1[["death"]] - by_hand) < 1e-12, d1[["death"]] > 0)   # the prose says 'ca tử vong thêm'

Bài 2 (trung bình). Ở ngưỡng của khóa học, giá xét nghiệm tối đa đổi thế nào khi tỷ lệ kháng là 10%, 20%, mức cơ sở (30.2%) và 50%? Tính hai cách: chỉ so với thường quy, và so với cả bốn chiến lược.

Chỉ so với thường quy, giá tối đa tăng từ 15.81 đến 80.77 triệu VNĐ. So với cả bốn chiến lược nó không đơn điệu: 15.81 / 32.47 / 27.41 / 14.36 triệu VNĐ, vì khi kháng nhiều, đối thủ mạnh nhất chuyển từ thường quy sang phủ CR cho mọi ca. Câu trả lời phụ thuộc đối thủ nào đang đứng đầu.

# Exercise 2: maximum test price (million VND) by CR prevalence at the course threshold:
# against usual care only, and against the best of all four strategies
pcr_values <- c(0.10, 0.20, p$p_cr, 0.50)
price_by_pcr <- data.frame(p_cr = pcr_values, vs_usual_only = NA_real_, vs_all = NA_real_, competitor = NA_character_)
for (k in seq_along(pcr_values)) {
  q <- modifyList(p, list(p_cr = pcr_values[k]))
  price_by_pcr$vs_usual_only[k] <- threshold_value(q, "c_test", lower = -1e8, upper = 1e9, outcome = inmb)
  price_by_pcr$vs_all[k] <- threshold_value(q, "c_test", lower = -1e8, upper = 1e9)
  price_by_pcr$competitor[k] <- STRATEGY_LABELS[[best_other(modifyList(q, list(c_test = 1.01 * price_by_pcr$vs_all[k])))]]
}
price_by_pcr[, c("vs_usual_only", "vs_all")] <- round(price_by_pcr[, c("vs_usual_only", "vs_all")] / 1e6, 2)
price_by_pcr
   p_cr vs_usual_only vs_all        competitor
1 0.100         15.81  15.81        Thường quy
2 0.200         32.47  32.47        Thường quy
3 0.302         49.30  27.41 Phủ CR cho mọi ca
4 0.500         80.77  14.36 Phủ CR cho mọi ca
# Against all strategies the maximum price is not monotone in prevalence: it peaks inside
stopifnot(which.max(price_by_pcr$vs_all) %in% 2:3, all(diff(price_by_pcr$vs_usual_only) > 0))
# The strongest rival moves from usual care (low resistance) to cover-all (high resistance)
stopifnot(price_by_pcr$competitor[1] == STRATEGY_LABELS[["usual"]],
          price_by_pcr$competitor[4] == STRATEGY_LABELS[["cover_all"]])

Bài 3 (khó). Tính tử vong khi điều trị sớm bằng cách sai, chia xác suất cho tỷ số chênh. Cách sai này ngầm giả định tỷ số chênh bằng bao nhiêu?

Cách sai cho 33.6% thay vì 42.0%, tức ngầm giả định tỷ số chênh 2.30 thay vì 1.60: hiệu quả bị phóng đại và ICER bị kéo xuống.

# Exercise 3: the wrong way to apply an odds ratio (dividing the probability)
m_wrong <- inp$m_cr_late / p$or_inapprop
stopifnot(m_wrong < inp$m_cr_early)   # the prose says the wrong method exaggerates the benefit
round(c(correct = inp$m_cr_early, wrong = m_wrong,
        or_implied_by_wrong = odds(inp$m_cr_late) / odds(m_wrong)), 4)
            correct               wrong or_implied_by_wrong 
             0.4202              0.3356              2.2958

Tài liệu tham khảo

# Slide Lời giảng chính Hình/Code
1 Câu hỏi và cây Hai giờ sáng ở ICU: colistin ngay, chờ kháng sinh đồ hay làm PCR? Nút quyết định, nút cơ hội, nút kết thúc; vẽ trước khi code Sơ đồ Mermaid
2 Từ y văn đến xác suất Tách tử vong theo kháng thuốc; số ngày ICU; kiểm tra trên đầu ra của cây Chunk l05-derive
3 Xét nghiệm đo gì? Phết trực tràng đo định cư, không đo tác nhân gây bệnh; cov_cr Bảng chiến lược
4 Tỷ số chênh Áp lên odds, không chia xác suất Bài tập 3
5 Bảng đường đi và rollback Tích xác suất nhánh; tính ngược = tổng theo đường đi Chunks l05-paths, l05-rollback
6 Kết quả Chi phí, QALY, ICER; biên so với đối thủ mạnh nhất Chunk l05-results
7 OWSA đúng câu hỏi So với phương án tốt nhất còn lại; loại khoảng thấp–cao; không tham số nào làm đổi kết luận figs/l05-tornado.png
8 Ngưỡng hai phía Tỷ lệ kháng thấp nhất và cao nhất; giá tối đa theo đối thủ figs/l05-threshold.png
9 Kịch bản A. baumannii và xét nghiệm chỉ bắt một nửa ca: phủ CR thắng; hiệu quả bằng nửa: vẫn giữ kết luận Chunk l05-scenarios
10 Kiểm thử hộp đen Hai nhánh giống hệt thì kết quả giống hệt Bài tập 1
11 Khi cây gãy Sự kiện lặp lại; 2^n đường đi; cầu nối sang Markov (M02) Chunk l05-limits