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")
source("code/model.R")
options(warnPartialMatchDollar = TRUE) # `$` that silently matches a shortened name now warns
p_sim <- read_params("params/amr_sim.csv")
ref <- read_params("params/reference_data.csv")
d <- read_params("params/teaching_demos.csv")
# The simulated "truth" (only possible because the data are simulated)
iso_true <- simulate_isolates(p_sim, seed = 2025)
# The raw laboratory export is written by Lesson 1. Raw data are read-only, so this lesson
# re-creates the file (same seeds as Lesson 1) only if it is missing, and never overwrites it.
raw_file <- "data/amr_isolates_raw.csv"
if (!file.exists(raw_file)) {
dir.create("data", showWarnings = FALSE)
write.csv(make_messy(iso_true, p_sim, seed = 2026), raw_file, row.names = FALSE, na = "")
}Bài 4: Data frame — nhập, làm sạch, tóm tắt và chuyển dạng dữ liệu
Từ tệp trích xuất kháng sinh đồ lộn xộn đến bảng tỷ lệ đề kháng có khoảng tin cậy
Vì sao bạn nên quan tâm?
Mỗi dòng trong tệp xuất của khoa vi sinh là tiền thật: một lần cấy máu bằng hệ thống tự động (325,200 VNĐ) cộng một kháng sinh đồ (213,800 VNĐ) đã bằng 2.2 ngày lương tối thiểu vùng I theo mức từ 1/1/2026 (Bộ Y tế 2024; Chính phủ 2025). Tệp xuất đó thường lộn xộn: tên vi khuẩn gõ nhiều kiểu, kết quả viết thường, và một người bệnh nhiễm chủng kháng được cấy đi cấy lại. Nếu không loại chủng trùng, tỷ lệ kháng carbapenem bị kéo lên và phác đồ kinh nghiệm của bệnh viện có thể chuyển sang colistin sớm hơn cần thiết. Trong VINARES 2016–2017, 79% chủng Acinetobacter kháng carbapenem, tính trên chủng đã loại trùng (Vu và c.s. 2021). Con số như vậy chỉ đáng tin khi mọi bước làm sạch và loại trùng được viết thành mã; bài này dạy các bước đó.
Mục tiêu học tập
Sau bài này, bạn có thể:
- Nhập một tệp CSV vào data frame với các đối số đúng cho giá trị thiếu và mã hóa tiếng Việt.
- Phát hiện và sửa bằng mã các lỗi nhập liệu điển hình: tên viết nhiều kiểu, mã kết quả không chuẩn, ngày dạng chữ, dòng trùng.
- Giữ lại chủng phân lập (isolate) đầu tiên của mỗi người bệnh và giải thích vì sao bước này làm thay đổi tỷ lệ đề kháng.
- Chuyển dữ liệu từ dạng rộng (wide format) sang dạng dài (long format) và ngược lại.
- Tính tỷ lệ đề kháng kèm khoảng tin cậy Wilson, kiểm tra kết quả với hàm có sẵn của R, và nói rõ khoảng tin cậy đó đo loại bất định nào.
Ý chính
- Data frame là bảng có các cột cùng độ dài; mỗi cột là một vectơ có kiểu riêng. Lấy dòng và cột bằng
df[dòng, cột]. - Không sửa dữ liệu thô. Mỗi bước làm sạch là một dòng mã, kèm một kiểm tra: tên lạ thì dừng, mã lạ thì dừng.
- Dùng bảng tra (lookup table) để chuẩn hóa tên. Đây là một vectơ có tên, dễ đọc và dễ bổ sung hơn một chuỗi
ifelse(). - Báo cáo tỷ lệ đề kháng dùng một chủng cho mỗi người bệnh, mỗi loài, mỗi kỳ (Hindler và Stelling 2007). Giữ chủng lặp lại sẽ làm lệch tỷ lệ về phía những người bệnh được cấy nhiều lần.
- Dữ liệu dạng dài (một dòng cho mỗi chủng × kháng sinh) là dạng thuận tiện nhất để tóm tắt và vẽ biểu đồ (Wickham 2014).
- Khoảng tin cậy Wald (\(\hat p \pm 1.96\,SE\)) hoạt động kém khi mẫu nhỏ hoặc tỷ lệ gần 0 hay 1. Hãy dùng khoảng Wilson (Wilson 1927; Agresti và Coull 1998; Brown, Cai, và DasGupta 2001). Khoảng này chỉ phản ánh sai số lấy mẫu nhị thức; khi gộp các bệnh viện và năm có khác biệt thật rồi muốn khái quát ra ngoài các ô đã quan sát, khoảng thực sự rộng hơn.
- Định nghĩa phải nói rõ: “đề kháng” chỉ tính R, hay “không nhạy cảm” tính cả I và R. Coi I là không nhạy cảm là quy ước của Clinical and Laboratory Standards Institute (CLSI) và của Magiorakos; theo European Committee on Antimicrobial Susceptibility Testing (EUCAST) hiện hành, I nghĩa là nhạy cảm khi tăng liều.
Nội dung
Data frame và nhập dữ liệu
read.csv() đọc tệp CSV thành data frame. Ba đối số quan trọng: na.strings cho biết ô nào là giá trị thiếu (ở đây là ô trống); encoding = "UTF-8" giữ đúng dấu tiếng Việt; và từ R 4.0, chuỗi ký tự không còn bị tự đổi thành factor. Ngay sau khi nhập, luôn xem dim(), str() và table() của các cột phân loại. Phần lớn lỗi dữ liệu lộ ra ở bước này.
Làm sạch bằng mã
Tệp trích xuất từ phòng xét nghiệm trong module này có bốn loại lỗi thường gặp trong dữ liệu thật: tên vi khuẩn viết nhiều kiểu (E.coli, Escherichia coli, A. baumanii), mã kết quả viết thường hoặc thừa khoảng trắng, ngày ở dạng chữ dd/mm/yyyy, và dòng nhập trùng. Nguyên tắc chung: chuẩn hóa trước (trimws(), tolower(), toupper()), tra bảng sau, và dừng nếu gặp giá trị không có trong bảng. Không bao giờ để một tên lạ lặng lẽ trở thành NA rồi biến mất khỏi mẫu số.
Dữ liệu ở đây là mô phỏng, nên ta có một lợi thế không bao giờ có với dữ liệu thật: so sánh kết quả làm sạch với dữ liệu gốc trước khi bị làm “bẩn”. Đây là cách tốt để kiểm tra một quy trình làm sạch trước khi dùng nó cho dữ liệu thật.
Chủng đầu tiên
Một người bệnh nhiễm khuẩn kéo dài thường được cấy nhiều lần, và các chủng lặp lại gần như giống hệt nhau. Nếu đếm tất cả, tỷ lệ đề kháng phản ánh số lần cấy chứ không phản ánh số người bệnh. Hướng dẫn của CLSI về bản tổng kết mức độ nhạy cảm kháng sinh tích lũy, bản M39-A2 mà Hindler và Stelling mô tả, coi việc loại chủng lặp lại của cùng một người bệnh là một vấn đề then chốt (Hindler và Stelling 2007). Mạng lưới giám sát ở Việt Nam cũng báo cáo trên chủng đã loại trùng: VINARES (Vu và c.s. 2021), và ACORN2 giữ chủng đầu tiên của mỗi loài theo loại bệnh phẩm trong mỗi đợt nhiễm khuẩn (Nguyen và c.s. 2026). Bài này dùng quy tắc đơn giản: chủng đầu tiên của mỗi người bệnh, mỗi loài, mỗi năm.
Dữ liệu rộng và dữ liệu dài
Tệp trích xuất có dạng rộng (wide): mỗi chủng một dòng, mỗi kháng sinh một cột. Để tóm tắt theo loài × kháng sinh, ta cần dạng dài (long): mỗi dòng là một cặp chủng × kháng sinh có kết quả. Dạng dài là dữ liệu “gọn” (tidy): mỗi biến một cột, mỗi quan sát một dòng (Wickham 2014). Hàm to_long() trong code/model.R chuyển dạng bằng base R; tidyr::pivot_longer() làm cùng việc đó, và sách R for Data Science hướng dẫn chi tiết cách làm theo tidyverse (Wickham, Çetinkaya-Rundel, và Grolemund không ngày). Bảng báo cáo cuối cùng lại có dạng rộng (loài × kháng sinh), tạo bằng tapply().
Tỷ lệ và khoảng tin cậy
Với \(x\) chủng đề kháng trên \(n\) chủng có kết quả, \(\hat p = x/n\). Khoảng Wald \(\hat p \pm z\sqrt{\hat p(1-\hat p)/n}\) có thể vượt ra ngoài [0, 1] và có độ rộng bằng 0 khi \(x = 0\). Khoảng Wilson (Wilson 1927) dời tâm về phía 0.5 và luôn nằm trong [0, 1]:
\[ \tilde p = \frac{\hat p + \dfrac{z^2}{2n}}{1 + \dfrac{z^2}{n}}, \qquad \text{nửa độ rộng} = \frac{z}{1 + z^2/n}\sqrt{\frac{\hat p(1-\hat p)}{n} + \frac{z^2}{4n^2}}, \]
với \(z = 1.96\) cho KTC 95%. Các so sánh hệ thống cho thấy khoảng Wilson đạt tỷ lệ bao phủ gần mức danh nghĩa tốt hơn khoảng Wald, còn khoảng Clopper-Pearson “chính xác” thường rộng hơn cần thiết (Agresti và Coull 1998; Brown, Cai, và DasGupta 2001). Trong R, prop.test(x, n, correct = FALSE) trả về đúng khoảng Wilson, nên ta dùng nó để kiểm tra hàm wilson_ci() tự viết.
Một điều cần nói rõ là bảng tóm tắt ước lượng cái gì. Bảng trong bài này là tỷ lệ gộp 5 năm (2021–2025) của 3 bệnh viện, không phải tỷ lệ của một quần thể đồng nhất. Khoảng Wilson giả định mọi chủng là một mẫu ngẫu nhiên độc lập từ cùng một quần thể, nên nó chỉ đo sai số lấy mẫu nhị thức. Khi các bệnh viện và các năm thật sự khác nhau (dữ liệu mô phỏng được thiết kế như vậy, giống dữ liệu thật), và ta muốn nói về mức đề kháng của các bệnh viện và năm nói chung chứ không chỉ riêng 15 ô bệnh viện × năm này, khoảng gộp quá hẹp. Bản tổng kết tích lũy theo CLSI cũng được lập cho từng cơ sở y tế và thường theo năm (Hindler và Stelling 2007), chứ không gộp nhiều cơ sở và nhiều năm như bảng ở đây. Thực hành bên dưới đo mức độ quá hẹp này bằng độ phân tán (dispersion) của mô hình quasi-binomial.
Thực hành
Nạp hàm, các bảng tham số và tệp thô do Bài 1 tạo ra. Nếu tệp chưa có, script tạo lại bằng đúng mã của Bài 1; nếu đã có, script chỉ đọc, không bao giờ ghi đè. Nhờ vậy bài này chạy độc lập mà không làm thay đổi dữ liệu thô.
Nhập tệp CSV: ô trống là giá trị thiếu, mã hóa UTF-8. Kiểm tra: cột date có kiểu chữ (chr), chưa phải ngày.
raw <- read.csv(raw_file, na.strings = "", encoding = "UTF-8")
dim(raw)[1] 4100 11str(raw)'data.frame': 4100 obs. of 11 variables:
$ isolate_id: chr "I00001" "I00002" "I00003" "I00004" ...
$ patient_id: chr "P00005" "P00073" "P00190" "P00218" ...
$ hospital : chr "A" "A" "A" "A" ...
$ ward : chr "non-ICU" "ICU" "non-ICU" "ICU" ...
$ date : chr "01/01/2021" "01/01/2021" "01/01/2021" "01/01/2021" ...
$ organism : chr "A. baumannii" "E. coli" "E. coli" "E. coli" ...
$ CRO : chr NA " r" "R" "R" ...
$ CIP : chr "R" "R" " r" "R" ...
$ GEN : chr "R" "R" "R" "R" ...
$ SXT : chr "R" "R" "R" "R" ...
$ MEM : chr "R" "S" "S" "R" ...Tệp trích xuất có 4,100 dòng và 11 cột.
Khám phá trước khi sửa. table() của cột tên vi khuẩn và cột meropenem cho thấy ngay các kiểu viết lạ; duplicated() đếm dòng trùng hoàn toàn. Kiểm tra: ghi lại từng loại lỗi trước khi viết mã sửa.
sort(table(raw$organism), decreasing = TRUE) # several spellings of the same species
E. coli A. baumannii K. pneumoniae
1937 980 872
Escherichia coli E.coli E. COLI
56 48 40
e. coli k. pneumoniae A. baumanii
38 29 25
A.baumannii Acinetobacter baumannii Klebsiella pneumoniae
25 18 18
K.pneumoniae
14 table(raw$MEM, useNA = "ifany") # " r", " s": lower case and spaces
i r s I R S <NA>
3 63 100 134 1356 2351 93 sum(duplicated(raw)) # rows entered twice[1] 20head(raw$date, 3) # text, not dates[1] "01/01/2021" "01/01/2021" "01/01/2021"Cột tên vi khuẩn có 13 cách viết khác nhau cho ba loài, và có 20 dòng trùng hoàn toàn.
Làm sạch từng bước. Mỗi bước có một stopifnot(): mọi tên phải nằm trong bảng tra, mọi kết quả phải là S, I, R hoặc NA, mọi ngày phải đọc được. Kiểm tra: sau bước 4 chỉ còn đúng ba loài.
clean <- raw
# 1. Organism: lower-case, trim spaces, then look up the standard name
clean$organism <- unname(ORGANISM_LOOKUP[tolower(trimws(clean$organism))])
stopifnot(!anyNA(clean$organism)) # every spelling was recognised
# 2. Results: trim spaces and upper-case; only S, I, R or NA may remain
for (a in names(ANTIBIOTICS)) clean[[a]] <- toupper(trimws(clean[[a]]))
stopifnot(all(unlist(clean[names(ANTIBIOTICS)]) %in% c("S", "I", "R", NA)))
# 3. Dates: dd/mm/yyyy text -> Date, then derive the year
clean$date <- as.Date(clean$date, format = "%d/%m/%Y")
stopifnot(!anyNA(clean$date))
clean$year <- as.integer(format(clean$date, "%Y"))
# 4. Exact duplicate rows
n_before <- nrow(clean)
clean <- clean[!duplicated(clean), ]
c(rows_before = n_before, rows_after = nrow(clean), removed = n_before - nrow(clean))rows_before rows_after removed
4100 4080 20 table(clean$organism)
A. baumannii E. coli K. pneumoniae
1044 2112 924 Đã loại 20 dòng trùng.
Các bước trên được đóng gói trong clean_isolates(). Ta kiểm tra hàm cho cùng kết quả, và vì dữ liệu là mô phỏng, kiểm tra luôn rằng dữ liệu sau làm sạch trùng khớp hoàn toàn với dữ liệu gốc. Kiểm tra: hai stopifnot() đều qua. Phép so sánh với dữ liệu gốc cũng xác nhận rằng tệp trong data/ đúng là tệp do mã tạo ra.
# The packaged function gives the same result; with SIMULATED data we can also check
# against the truth, which is never possible with real data. This also confirms that the
# file in data/ is the export the code generates (nobody edited it).
clean_fn <- clean_isolates(raw)
stopifnot(isTRUE(all.equal(clean_fn, iso_true)))
stopifnot(nrow(clean) == nrow(clean_fn))Lấy một phần data frame theo điều kiện trên dòng và theo tên cột: các chủng A. baumannii từ ICU. Kiểm tra: chỉ còn bốn cột đã chọn.
# Rows by a condition, columns by name
icu_aba <- clean_fn[clean_fn$ward == "ICU" & clean_fn$organism == "A. baumannii",
c("isolate_id", "hospital", "year", "MEM")]
head(icu_aba, 4) isolate_id hospital year MEM
22 I00022 A 2021 R
24 I00024 A 2021 R
26 I00026 A 2021 R
45 I00045 A 2021 Rnrow(icu_aba)[1] 232Giữ chủng đầu tiên của mỗi người bệnh × loài × năm. Kiểm tra: không còn cặp người bệnh-loài-năm nào bị lặp.
iso <- clean_fn
first <- first_isolates(iso)
c(isolates = nrow(iso), first_isolates = nrow(first), patients = length(unique(iso$patient_id))) isolates first_isolates patients
4080 3094 3094 stopifnot(!anyDuplicated(first[, c("patient_id", "organism", "year")]))
# In these simulated data each patient has one species in one year, so first isolates = patients
stopifnot(nrow(first) == length(unique(iso$patient_id)))
table(first$organism)
A. baumannii E. coli K. pneumoniae
701 1713 680 Từ 4,080 chủng còn 3,094 chủng đầu tiên, gồm 701 chủng A. baumannii, 1,713 chủng E. coli và 680 chủng K. pneumoniae. Ở đây số chủng đầu tiên bằng số người bệnh, vì trong dữ liệu mô phỏng mỗi người bệnh chỉ có một loài trong một năm; stopifnot() trong đoạn mã đã kiểm tra điều đó.
Chuyển sang dạng dài. Kiểm tra: số dòng của bảng dài bằng số ô có kết quả (không phải NA) trong bảng rộng.
long <- to_long(first)
dim(long)[1] 14321 9head(long[, c("isolate_id", "organism", "antibiotic", "result")], 6) isolate_id organism antibiotic result
1 I00002 E. coli CRO R
2 I00003 E. coli CRO R
3 I00004 E. coli CRO R
4 I00005 K. pneumoniae CRO S
5 I00006 E. coli CRO R
6 I00007 E. coli CRO R# Check: one long row per non-missing result in the wide table
stopifnot(nrow(long) == sum(!is.na(first[, names(ANTIBIOTICS)])))Bảng dài có 14,321 dòng, mỗi dòng là một kết quả kháng sinh đồ.
Tóm tắt: số chủng đề kháng, số chủng có kết quả, tỷ lệ và KTC 95% Wilson cho mỗi loài × kháng sinh. Kiểm tra: stopifnot() xác nhận lo ≤ est ≤ hi ở mọi dòng, và không có dòng A. baumannii × ceftriaxone vì kháng sinh này không được thử nghiệm. Đây là tỷ lệ gộp mọi năm và mọi bệnh viện.
res <- resistance_summary(long, by = c("organism", "antibiotic"))
res_print <- res
res_print[, c("est", "lo", "hi")] <- round(res_print[, c("est", "lo", "hi")], 3)
res_print organism antibiotic n_r n est lo hi
1 A. baumannii CIP 562 683 0.823 0.792 0.850
2 E. coli CIP 1114 1659 0.671 0.649 0.694
3 K. pneumoniae CIP 312 666 0.468 0.431 0.506
4 E. coli CRO 1132 1644 0.689 0.666 0.710
5 K. pneumoniae CRO 352 655 0.537 0.499 0.575
6 A. baumannii GEN 499 682 0.732 0.697 0.764
7 E. coli GEN 803 1670 0.481 0.457 0.505
8 K. pneumoniae GEN 293 655 0.447 0.410 0.486
9 A. baumannii MEM 448 688 0.651 0.615 0.686
10 E. coli MEM 171 1660 0.103 0.089 0.119
11 K. pneumoniae MEM 228 666 0.342 0.307 0.379
12 A. baumannii SXT 658 682 0.965 0.948 0.976
13 E. coli SXT 1176 1652 0.712 0.690 0.733
14 K. pneumoniae SXT 355 659 0.539 0.501 0.576stopifnot(all(res$lo <= res$est & res$est <= res$hi)) # the estimate lies inside its own interval
range(res$n) # isolates behind each proportion[1] 655 1670Dữ liệu được mô phỏng quanh các mốc công bố của Việt Nam (các dòng target_* trong params/amr_sim.csv): VINARES 2016–2017 cho ceftriaxone, ciprofloxacin, gentamicin và trimethoprim-sulfamethoxazole (Vu và c.s. 2021), khảo sát bệnh viện toàn quốc cho carbapenem (Larsson và c.s. 2025). Đoạn mã sau so bảng vừa tính với các mốc đó. Cột z là chênh lệch chia cho sai số chuẩn nhị thức tại mốc. Kiểm tra: mọi giá trị z nằm trong khoảng ±4, tức chênh lệch còn trong tầm dao động lấy mẫu.
# The data are simulated AROUND published Vietnamese figures (target_* rows in params/amr_sim.csv).
# Compare the pooled proportions with their targets; z = difference / binomial SE at the target.
target <- resistance_matrix(p_sim)
res$target <- target[cbind(match(res$organism, ORGANISMS), match(res$antibiotic, colnames(target)))]
res$z <- (res$est - res$target) / sqrt(res$target * (1 - res$target) / res$n)
anchor_tab <- res[, c("organism", "antibiotic", "n", "est", "target", "z")]
anchor_tab[, c("est", "target", "z")] <- round(anchor_tab[, c("est", "target", "z")], 3)
anchor_tab organism antibiotic n est target z
1 A. baumannii CIP 683 0.823 0.816 0.461
2 E. coli CIP 1659 0.671 0.670 0.129
3 K. pneumoniae CIP 666 0.468 0.440 1.480
4 E. coli CRO 1644 0.689 0.717 -2.560
5 K. pneumoniae CRO 655 0.537 0.556 -0.958
6 A. baumannii GEN 682 0.732 0.738 -0.376
7 E. coli GEN 1670 0.481 0.477 0.314
8 K. pneumoniae GEN 655 0.447 0.465 -0.907
9 A. baumannii MEM 688 0.651 0.640 0.610
10 E. coli MEM 1660 0.103 0.110 -0.910
11 K. pneumoniae MEM 666 0.342 0.370 -1.478
12 A. baumannii SXT 682 0.965 0.957 1.005
13 E. coli SXT 1652 0.712 0.727 -1.381
14 K. pneumoniae SXT 659 0.539 0.524 0.755stopifnot(!anyNA(res$target), all(abs(res$z) < 4)) # within sampling error of the anchorsChênh lệch lớn nhất là 2.6 sai số chuẩn, ở cặp E. coli × CRO. Đây là lợi thế của dữ liệu mô phỏng: ta biết “sự thật” mà nó được sinh ra quanh đó, nên kiểm tra được cả quy trình tóm tắt.
Mỗi tỷ lệ trong bảng dựa trên từ 655 đến 1,670 chủng. Khoảng Wilson trong bảng trên chưa tính đến khác biệt giữa bệnh viện và giữa các năm. Để đo, ta chia dữ liệu meropenem thành các ô bệnh viện × năm, rồi ước lượng độ phân tán bằng glm() với họ quasibinomial (hàm có sẵn của base R; ở đây ta chỉ lấy con số dispersion, không dùng hệ số). Độ phân tán bằng 1 nghĩa là chỉ có sai số nhị thức; lớn hơn 1 nghĩa là các ô khác nhau nhiều hơn mức ngẫu nhiên cho phép. Kiểm tra: độ phân tán lớn hơn 1 ở cả ba loài.
# The Wilson CI treats all pooled isolates as ONE random sample. Do hospitals and years differ?
# Meropenem: resistance in each hospital x year cell, then a quasi-binomial model with no
# covariates. Its dispersion is 1 if only binomial sampling error is present, and above 1 if
# the cells differ more than chance allows.
mem_cells <- resistance_summary(long[long$antibiotic == "MEM", ],
by = c("organism", "hospital", "year"))
disp <- sapply(unique(mem_cells$organism), function(o) {
cells <- mem_cells[mem_cells$organism == o, ]
fit <- glm(cbind(n_r, n - n_r) ~ 1, family = quasibinomial, data = cells)
summary(fit)$dispersion
})
n_cells <- table(mem_cells$organism)
round(rbind(cells = n_cells[names(disp)], dispersion = disp, ci_wider_by = sqrt(disp)), 2) A. baumannii E. coli K. pneumoniae
cells 15.00 15.00 15.00
dispersion 2.78 2.66 2.57
ci_wider_by 1.67 1.63 1.60stopifnot(all(disp > 1))Với meropenem, độ phân tán nằm từ 2.6 đến 2.8, cao nhất ở A. baumannii (2.8). Khoảng tin cậy thực sự rộng hơn khoảng Wilson gộp xấp xỉ \(\sqrt{\text{độ phân tán}}\) lần, tức từ 1.6 đến 1.7 lần. Vì vậy hãy đọc khoảng Wilson của bảng gộp như một cận dưới của bất định. Muốn nói về một bệnh viện hay một năm, tóm tắt riêng từng nhóm đó (xem Bài 5), hoặc dùng mô hình có hiệu ứng nhóm.
So sánh với các hàm có sẵn, dùng ví dụ K. pneumoniae × meropenem. Kiểm tra: khoảng Wilson tự viết trùng với prop.test() đến 10 chữ số thập phân; khoảng Clopper-Pearson của binom.test() hơi rộng hơn. Hai dòng cuối cho thấy khoảng Wald hỏng ở mẫu nhỏ.
# Wilson interval = prop.test() without continuity correction; exact = binom.test()
k <- which(res$organism == "K. pneumoniae" & res$antibiotic == "MEM")
x <- res$n_r[k]; n <- res$n[k]
ci_ours <- c(res$lo[k], res$hi[k])
ci_prop <- as.vector(prop.test(x, n, correct = FALSE)$conf.int)
ci_exact <- as.vector(binom.test(x, n)$conf.int)
stopifnot(max(abs(ci_ours - ci_prop)) < 1e-10)
round(rbind(wilson = ci_ours, prop_test = ci_prop, clopper_pearson = ci_exact), 4) [,1] [,2]
wilson 0.3073 0.3792
prop_test 0.3073 0.3792
clopper_pearson 0.3063 0.3798# Why not Wald? With a few resistant isolates out of a small sample its lower limit is below zero
p_small <- d$wald_demo_x / d$wald_demo_n
wald_lo_small <- p_small - qnorm(0.975) * sqrt(p_small * (1 - p_small) / d$wald_demo_n)
wilson_lo_small <- wilson_ci(d$wald_demo_x, d$wald_demo_n)$lo
c(wald_lower = wald_lo_small, wilson_lower = wilson_lo_small) wald_lower wilson_lower
-0.03869417 0.03736126 Tỷ lệ K. pneumoniae đề kháng meropenem (%) là 34.2 (KTC 95%: 30.7–37.9), tính trên 666 chủng đầu tiên. Với 2 chủng đề kháng trên 15, cận dưới của khoảng Wald là -0.039, một tỷ lệ âm vô nghĩa; cận dưới của khoảng Wilson là 0.037.
Loại chủng lặp lại có đáng không? So sánh tỷ lệ đề kháng meropenem tính trên mọi chủng với tỷ lệ tính trên chủng đầu tiên. Kiểm tra: hai bảng xếp loài theo cùng thứ tự trước khi ghép.
# Meropenem resistance with all isolates versus first isolates only
mem_all <- resistance_summary(to_long(iso), by = c("organism", "antibiotic"))
mem_all <- mem_all[mem_all$antibiotic == "MEM", ]
mem_first <- res[res$antibiotic == "MEM", ]
dedup_tab <- data.frame(organism = mem_first$organism,
all_isolates = round(mem_all$est, 3),
first_isolates = round(mem_first$est, 3))
dedup_tab organism all_isolates first_isolates
1 A. baumannii 0.724 0.651
2 E. coli 0.136 0.103
3 K. pneumoniae 0.426 0.342stopifnot(identical(mem_all$organism, mem_first$organism))
# The lesson (and its opening box) says repeats pull the meropenem proportion UP for every species
stopifnot(all(mem_all$est > mem_first$est))
gap_kpn <- mem_all$est[mem_all$organism == "K. pneumoniae"] -
mem_first$est[mem_first$organism == "K. pneumoniae"]Với K. pneumoniae, tính trên mọi chủng cho 42.6%, còn tính trên chủng đầu tiên cho 34.2%: chênh 8.3 điểm phần trăm. Chênh lệch này có trong dữ liệu vì ta đã mô phỏng để người bệnh có chủng kháng meropenem được cấy lại thường hơn (xác suất cấy lại 35% so với 10% mỗi lần, giá trị minh họa).
Chuyển bảng tóm tắt dài thành bảng báo cáo rộng (loài × kháng sinh). Kiểm tra: ô A. baumannii × CRO là NA, không phải 0.
# Long summary -> wide antibiogram table (organism x antibiotic); NA = not tested
abg <- tapply(res$est, list(res$organism, res$antibiotic), sum)
abg <- abg[, names(ANTIBIOTICS)]
round(100 * abg) CRO CIP GEN SXT MEM
A. baumannii NA 82 73 96 65
E. coli 69 67 48 71 10
K. pneumoniae 54 47 45 54 34Cùng việc chuyển dạng và tóm tắt, viết bằng tidyr và dplyr. Nhiều người thấy cách viết này dễ đọc hơn; điều quan trọng là hai cách cho cùng một kết quả. Kiểm tra: mọi số đếm khớp nhau.
# The same reshaping and summary with tidyr/dplyr; results must match exactly
long_tv <- tidyr::pivot_longer(first, cols = c(CRO, CIP, GEN, SXT, MEM),
names_to = "antibiotic", values_to = "result",
values_drop_na = TRUE)
sum_tv <- long_tv |>
dplyr::group_by(organism, antibiotic) |>
dplyr::summarise(n_r = sum(result == "R"), n = dplyr::n(), .groups = "drop")
check <- merge(res[, c("organism", "antibiotic", "n_r", "n")], sum_tv,
by = c("organism", "antibiotic"), suffixes = c("_base", "_tv"))
stopifnot(nrow(check) == nrow(res), all(check$n_r_base == check$n_r_tv),
all(check$n_base == check$n_tv))
head(check, 3) organism antibiotic n_r_base n_base n_r_tv n_tv
1 A. baumannii CIP 562 683 562 683
2 A. baumannii GEN 499 682 499 682
3 A. baumannii MEM 448 688 448 688
Mẹo thực hành
- Đọc dữ liệu một lần ở đầu script, lưu thành đối tượng
raw, và không bao giờ ghi đè lên nó. Mọi bước làm sạch tạo đối tượng mới. - Viết bảng tra tên vi khuẩn thành một vectơ có tên trong
code/model.R. Khi phòng xét nghiệm có kiểu viết mới,clean_isolates()sẽ dừng và nêu đúng tên đó; bạn chỉ cần thêm một dòng vào bảng tra. - Khi so sánh hai bảng tóm tắt, dùng
merge()theo khóa (organism,antibiotic) thay vì giả định hai bảng có cùng thứ tự dòng. - Trong
$, luôn gõ đủ tên cột, và tên phần tử nếu là danh sách (Bài 2). R cho phép gõ tắt: nếu tên tắt khớp đúng một cột, R trả về cột đó; nếu khớp nhiều cột, R trả vềNULLmà không báo lỗi. Cảnh báowarnPartialMatchDollarbắt được lỗi này.
Lỗi thường gặp
- Tính tỷ lệ trên mọi chủng. Tỷ lệ bị kéo về phía những người bệnh được cấy nhiều lần. Luôn nói rõ quy tắc loại trùng đã dùng.
- Đổi
NAthành"S". Chủng chưa thử nghiệm không phải chủng nhạy cảm.NAphải nằm ngoài cả tử số lẫn mẫu số. - Dùng sai khoảng tin cậy. Khoảng Wald cho mẫu nhỏ có cận dưới âm như ví dụ ở trên; hãy dùng
wilson_ci()hoặcprop.test(correct = FALSE). Và đừng coi khoảng Wilson của bảng gộp là toàn bộ bất định: nó không tính khác biệt giữa bệnh viện và giữa các năm. - Kiểm tra không bao giờ thất bại.
all(logical(0))làTRUE, nênstopifnot()trên một vectơ rỗng luôn qua. Điều này dễ xảy ra khi$trả vềNULLvì tên viết tắt không khớp cột nào: ví dụcheck$ncó thể khớp nhầm hoặc không khớp. Hãy gõ đủ tên (check$n_base) và thử kiểm tra trên một ví dụ đã biết là sai.
Bài tập
Bài tập 1 (dễ). Với K. pneumoniae và gentamicin, so sánh tỷ lệ “đề kháng” (chỉ R) với tỷ lệ “không nhạy cảm” (I hoặc R), kèm KTC 95%.
Lời giải
Đối số nonsusceptible = TRUE của resistance_summary() đổi tử số từ R sang I + R.
# Exercise 1: K. pneumoniae and gentamicin: resistant (R) vs non-susceptible (I + R)
r_only <- resistance_summary(long, nonsusceptible = FALSE)
ns <- resistance_summary(long, nonsusceptible = TRUE)
k1 <- r_only$organism == "K. pneumoniae" & r_only$antibiotic == "GEN"
round(rbind(R = unlist(r_only[k1, c("est", "lo", "hi")]),
I_or_R = unlist(ns[k1, c("est", "lo", "hi")])), 3) est lo hi
R 0.447 0.410 0.486
I_or_R 0.476 0.438 0.515Bài tập 2 (vừa). Tính tỷ lệ đề kháng meropenem theo loài và theo khoa (ICU so với khoa khác), dùng chủng đầu tiên. Ở loài nào KTC của hai khoa không chồng lên nhau?
Lời giải
Lọc bảng dài theo kháng sinh rồi tóm tắt theo organism và ward. Hai KTC không chồng lên nhau là dấu hiệu khác biệt rõ; hai KTC chồng lên nhau chưa đủ để kết luận là không khác biệt.
# Exercise 2: meropenem resistance by ward (ICU vs other wards), first isolates
by_ward <- resistance_summary(long[long$antibiotic == "MEM", ], by = c("organism", "ward"))
by_ward[, c("est", "lo", "hi")] <- round(by_ward[, c("est", "lo", "hi")], 3)
by_ward organism ward n_r n est lo hi
1 A. baumannii ICU 108 146 0.740 0.663 0.804
2 E. coli ICU 47 312 0.151 0.115 0.195
3 K. pneumoniae ICU 57 127 0.449 0.365 0.536
4 A. baumannii non-ICU 340 542 0.627 0.586 0.667
5 E. coli non-ICU 124 1348 0.092 0.078 0.109
6 K. pneumoniae non-ICU 171 539 0.317 0.279 0.358Bài tập 3 (khó hơn, dữ liệu thật). Báo cáo VINARES 2016–2017 công bố số chủng Acinetobacter spp. và Pseudomonas aeruginosa kháng carbapenem trên số chủng có kết quả (Vu và c.s. 2021). Các con số này nằm trong params/reference_data.csv. Tính KTC 95% Wilson cho hai tỷ lệ. Vì sao KTC hẹp hơn nhiều so với KTC của bệnh viện mô phỏng?
Lời giải
# Exercise 3: real data. Carbapenem resistance in VINARES 2016-2017 (vu2021antimicrobial)
vinares <- rbind(
acinetobacter = wilson_ci(ref$vinares_acb_carb_r, ref$vinares_acb_carb_n),
pseudomonas = wilson_ci(ref$vinares_pae_carb_r, ref$vinares_pae_carb_n)
)
round(vinares, 3) est lo hi
acinetobacter 0.788 0.775 0.801
pseudomonas 0.448 0.432 0.465Tỷ lệ Acinetobacter spp. kháng carbapenem trong VINARES là 78.8 (KTC 95%: 77.5–80.1) (%). KTC hẹp vì mẫu số lớn: 3,622 chủng, so với 655 đến 1,670 chủng ở mỗi loài × kháng sinh trong dữ liệu mô phỏng (KTC của K. pneumoniae × meropenem rộng 7.2 điểm phần trăm, so với 2.7 của VINARES). KTC chỉ phản ánh sai số chọn mẫu ngẫu nhiên; nó không phản ánh sai lệch do cách chọn bệnh phẩm hay do khác biệt giữa các bệnh viện, mà chính các tác giả VINARES đã lưu ý.
Tài liệu tham khảo
Kịch bản video (~15 phút)
| # | Slide | Lời giảng chính | Hình/Code |
|---|---|---|---|
| 1 | Data frame | Cột là vectơ; df[dòng, cột]; str() |
l04-import |
| 2 | Khám phá trước khi sửa | table(), duplicated(); bốn loại lỗi |
l04-inspect |
| 3 | Làm sạch bằng mã | Chuẩn hóa → tra bảng → dừng nếu lạ | l04-clean-steps |
| 4 | Kiểm tra quy trình làm sạch | Dữ liệu mô phỏng cho phép so với “sự thật” | l04-clean-check |
| 5 | Chủng đầu tiên và mốc công bố | Vì sao loại trùng; CLSI M39-A2; VINARES, ACORN2; dữ liệu mô phỏng quanh mốc công bố, cột z |
l04-first-isolate, l04-anchor-check |
| 6 | Rộng và dài | Dữ liệu gọn; to_long(); pivot_longer() |
l04-long |
| 7 | Khoảng Wilson | Wald hỏng ở mẫu nhỏ; công thức; kiểm tra với prop.test(); ước lượng gì (gộp mọi năm và bệnh viện) và độ phân tán |
l04-ci-check, l04-dispersion |
| 8 | Ảnh hưởng của loại trùng | Mọi chủng so với chủng đầu tiên | l04-dedup-effect |
| 9 | Bảng báo cáo | Dài → rộng bằng tapply(); NA khác 0 |
l04-wide |
| 10 | Cách viết tidyverse | Cùng kết quả, kiểm tra bằng merge() |
l04-tidyverse |
| 11 | Dữ liệu thật | KTC cho số liệu VINARES; giới hạn của KTC | l04-ex3 |