RMarkdown_DATA_AI_POLYP
Đào Anh Dũng
2026-09-30
PHẦN MỞ ĐẦU. THIẾT LẬP MÔI TRƯỜNG LÀM VIỆC
0.1. Thiết lập thư mục làm việc
Thư mục làm việc (working directory) là nơi R tìm file dữ liệu và ghi file kết quả. Theo yêu cầu, dùng lệnh setwd() với đường dẫn tuyệt đối tới thư mục Data.
Lưu ý kỹ thuật quan trọng: khi bấm Knit, knitr đặt lại thư mục làm việc về thư mục chứa file .Rmd sau mỗi khối lệnh, nên setwd() chỉ có tác dụng trong khối chứa nó. Vì vậy phải thêm knitr::opts_knit$set(root.dir = …) với cùng đường dẫn – hai dòng này đi cùng nhau.
# Đường dẫn tuyệt đối tới thư mục Data (thư mục gốc làm việc)
thu_muc_goc <- "/Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Data"
if (!dir.exists(thu_muc_goc)) stop("Không tìm thấy thư mục: ", thu_muc_goc,
"\nKiểm tra OneDrive đã đồng bộ và đường dẫn còn đúng.")
setwd(thu_muc_goc) # khi chạy từng khối trong RStudio
knitr::opts_knit$set(root.dir = thu_muc_goc) # khi bấm Knit
# Thư mục lưu kết quả: …/AI_Colonoscopy_STARD/Output/ket_qua (tính từ thư mục Data)
thu_muc_kq <- normalizePath(file.path(thu_muc_goc, "..", "Output", "ket_qua"), mustWork = FALSE)
dir.create(thu_muc_kq, recursive = TRUE, showWarnings = FALSE)
file_du_lieu <- file.path(thu_muc_goc, "DATA_AI_POLYP.xlsx")
sheet_du_lieu <- "C"
if (!file.exists(file_du_lieu)) stop("Không tìm thấy file dữ liệu: ", file_du_lieu)
cat("Thư mục làm việc :", getwd(), "\nFile dữ liệu :", file_du_lieu,
"\nThư mục kết quả :", thu_muc_kq, "\n")## Thư mục làm việc : /Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Data
## File dữ liệu : /Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Data/DATA_AI_POLYP.xlsx
## Thư mục kết quả : /Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Output/ket_qua# Đường dẫn tuyệt đối tới thư mục Data (thư mục gốc làm việc)
thu_muc_goc <- "/Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Data"
if (!dir.exists(thu_muc_goc)) stop("Không tìm thấy thư mục: ", thu_muc_goc,
"\nKiểm tra OneDrive đã đồng bộ và đường dẫn còn đúng.")
setwd(thu_muc_goc) # khi chạy từng khối trong RStudio
knitr::opts_knit$set(root.dir = thu_muc_goc) # khi bấm Knit
# Thư mục lưu kết quả: …/AI_Colonoscopy_STARD/Output/ket_qua (tính từ thư mục Data)
thu_muc_kq <- normalizePath(file.path(thu_muc_goc, "..", "Output", "ket_qua"), mustWork = FALSE)
dir.create(thu_muc_kq, recursive = TRUE, showWarnings = FALSE)
file_du_lieu <- file.path(thu_muc_goc, "DATA_AI_POLYP.xlsx")
sheet_du_lieu <- "C"
if (!file.exists(file_du_lieu)) stop("Không tìm thấy file dữ liệu: ", file_du_lieu)
cat("Thư mục làm việc :", getwd(), "\nFile dữ liệu :", file_du_lieu,
"\nThư mục kết quả :", thu_muc_kq, "\n")## Thư mục làm việc : /Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Data
## File dữ liệu : /Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Data/DATA_AI_POLYP.xlsx
## Thư mục kết quả : /Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Output/ket_quasetwd("/Users/daoanhdung/Library/CloudStorage/OneDrive-Personal/NCKH/AI_ColonoScopy_Polyp/AI_Colonoscopy_STARD/Data")0.2. Cài đặt packages
danh_sach_goi <- data.frame(
`Gói` = c("rmarkdown", "knitr", "readxl", "dplyr", "tidyr", "stringi",
"compareGroups", "table1", "flextable", "officer", "writexl",
"geepack", "DTComPair",
"ggplot2", "ggpubr", "ggsignif", "patchwork", "scales", "ragg", "systemfonts"),
`Nhóm việc` = c("Báo cáo", "Báo cáo", "Nhập dữ liệu", "Xử lý dữ liệu", "Xử lý dữ liệu", "Xử lý dữ liệu",
"Bảng mô tả", "Bảng mô tả", "Bảng Word", "Bảng Word", "Xuất Excel",
"Thống kê", "Thống kê",
"Biểu đồ", "Biểu đồ", "Biểu đồ", "Biểu đồ", "Biểu đồ", "Xuất hình", "Xuất hình"),
`Công dụng trong phân tích này` = c(
"Chuyển file .Rmd thành báo cáo HTML/Word (nút Knit)",
"Chạy các khối lệnh R và chèn kết quả vào báo cáo",
"Đọc file Excel DATA_AI_POLYP.xlsx",
"Lọc, tạo biến mới, gom nhóm, tổng hợp dữ liệu (mutate, filter, group_by, summarise)",
"Chuyển dữ liệu dạng rộng ↔ dài",
"Chuẩn hoá Unicode tiếng Việt (NFC) – tránh lỗi chữ có dấu từ Excel trên macOS",
"Bảng đặc điểm mẫu có p-value tự động (t-test/Mann-Whitney, χ²/Fisher); xuất Word, Excel",
"Bảng mô tả dạng HTML đẹp, có cột Chung (overall) và cột p tự viết",
"Tạo bảng định dạng Word: font Times New Roman 11, kẻ khung, tô màu tiêu đề",
"Thao tác file Word (flextable dùng kèm)",
"Ghi các bảng kết quả ra file Excel",
"Mô hình GEE – tính p hiệu chỉnh cụm theo bệnh nhân (1 bệnh nhân có nhiều polyp)",
"So sánh 2 test chẩn đoán bắt cặp: McNemar chính xác (Se, Sp), Leisenring (PPV, NPV), Gu–Pepe (LR)",
"Vẽ biểu đồ (chuẩn trình bày của GIE, Clinical Endoscopy)",
"Giao diện biểu đồ chuẩn tạp chí (theme_pubr) và gắn ngoặc p-value lên biểu đồ",
"Vẽ ngoặc so sánh và p-value giữa hai cột (được ggpubr sử dụng)",
"Ghép nhiều biểu đồ thành một hình (panel A, B)",
"Định dạng trục phần trăm",
"Xuất PNG 300 dpi và TIFF 600 dpi nén LZW, hiển thị đúng tiếng Việt",
"Tìm font Arial trên máy để dùng cho chữ trong hình"),
check.names = FALSE)
knitr::kable(danh_sach_goi, caption = "Các gói R sử dụng")| Gói | Nhóm việc | Công dụng trong phân tích này |
|---|---|---|
| rmarkdown | Báo cáo | Chuyển file .Rmd thành báo cáo HTML/Word (nút Knit) |
| knitr | Báo cáo | Chạy các khối lệnh R và chèn kết quả vào báo cáo |
| readxl | Nhập dữ liệu | Đọc file Excel DATA_AI_POLYP.xlsx |
| dplyr | Xử lý dữ liệu | Lọc, tạo biến mới, gom nhóm, tổng hợp dữ liệu (mutate, filter, group_by, summarise) |
| tidyr | Xử lý dữ liệu | Chuyển dữ liệu dạng rộng ↔︎ dài |
| stringi | Xử lý dữ liệu | Chuẩn hoá Unicode tiếng Việt (NFC) – tránh lỗi chữ có dấu từ Excel trên macOS |
| compareGroups | Bảng mô tả | Bảng đặc điểm mẫu có p-value tự động (t-test/Mann-Whitney, χ²/Fisher); xuất Word, Excel |
| table1 | Bảng mô tả | Bảng mô tả dạng HTML đẹp, có cột Chung (overall) và cột p tự viết |
| flextable | Bảng Word | Tạo bảng định dạng Word: font Times New Roman 11, kẻ khung, tô màu tiêu đề |
| officer | Bảng Word | Thao tác file Word (flextable dùng kèm) |
| writexl | Xuất Excel | Ghi các bảng kết quả ra file Excel |
| geepack | Thống kê | Mô hình GEE – tính p hiệu chỉnh cụm theo bệnh nhân (1 bệnh nhân có nhiều polyp) |
| DTComPair | Thống kê | So sánh 2 test chẩn đoán bắt cặp: McNemar chính xác (Se, Sp), Leisenring (PPV, NPV), Gu–Pepe (LR) |
| ggplot2 | Biểu đồ | Vẽ biểu đồ (chuẩn trình bày của GIE, Clinical Endoscopy) |
| ggpubr | Biểu đồ | Giao diện biểu đồ chuẩn tạp chí (theme_pubr) và gắn ngoặc p-value lên biểu đồ |
| ggsignif | Biểu đồ | Vẽ ngoặc so sánh và p-value giữa hai cột (được ggpubr sử dụng) |
| patchwork | Biểu đồ | Ghép nhiều biểu đồ thành một hình (panel A, B) |
| scales | Biểu đồ | Định dạng trục phần trăm |
| ragg | Xuất hình | Xuất PNG 300 dpi và TIFF 600 dpi nén LZW, hiển thị đúng tiếng Việt |
| systemfonts | Xuất hình | Tìm font Arial trên máy để dùng cho chữ trong hình |
goi_can_thiet <- c("rmarkdown", "knitr", "readxl", "dplyr", "tidyr", "stringi",
"compareGroups", "table1", "flextable", "officer", "writexl",
"geepack", "DTComPair",
"ggplot2", "ggpubr", "ggsignif", "patchwork", "scales", "ragg", "systemfonts")
goi_thieu <- goi_can_thiet[!goi_can_thiet %in% rownames(installed.packages())]
if (length(goi_thieu) > 0) {
cat("Đang cài các gói còn thiếu:", paste(goi_thieu, collapse = ", "), "\n")
install.packages(goi_thieu, repos = "https://cloud.r-project.org", dependencies = TRUE)
} else {
cat("Tất cả", length(goi_can_thiet), "gói đã được cài đặt – bỏ qua bước cài.\n")
}## Tất cả 20 gói đã được cài đặt – bỏ qua bước cài.0.3. Gọi thư viện (library)
suppressPackageStartupMessages({
library(readxl) # đọc Excel
library(dplyr) # xử lý dữ liệu
library(tidyr) # định dạng dữ liệu dài/rộng
library(compareGroups) # bảng đặc điểm mẫu + p-value
library(table1) # bảng mô tả HTML
library(flextable) # bảng Word
library(officer) # hỗ trợ flextable
library(geepack) # GEE – hiệu chỉnh cụm theo bệnh nhân
library(DTComPair) # so sánh 2 test chẩn đoán bắt cặp
library(ggplot2) # biểu đồ
library(ggpubr) # theme tạp chí + p-value trên biểu đồ
library(ggsignif) # ngoặc so sánh
library(patchwork) # ghép biểu đồ
library(scales) # định dạng %
})
knitr::opts_chunk$set(dev = "ragg_png", dpi = 150) # hiển thị hình trong báo cáo bằng ragg
cat("R", as.character(getRversion()), "| compareGroups", as.character(packageVersion("compareGroups")),
"| table1", as.character(packageVersion("table1")), "| DTComPair", as.character(packageVersion("DTComPair")),
"| ggplot2", as.character(packageVersion("ggplot2")), "\n")## R 4.6.1 | compareGroups 4.10.4 | table1 1.5.1 | DTComPair 1.2.6 | ggplot2 4.0.30.4. Tham số phân tích và các hàm định dạng dùng chung
B = 2000 lần bootstrap cụm theo bệnh nhân, hạt giống (seed) = 2026. Hai tham số này giống phiên phân tích trước, nên khoảng tin cậy (KTC) 95% sẽ trùng khớp.
Màu Okabe–Ito: an toàn cho người mù màu và vẫn phân biệt được khi in đen trắng. Font trong hình: Arial (không chân) theo thông lệ tạp chí; font trong bảng Word: Times New Roman 11.
Hàm dinh_dang_bang() định dạng mọi bảng giống nhau; hàm luu_hinh() lưu mỗi hình thành PNG 300 dpi + TIFF 600 dpi (LZW).
B <- 2000 # số lần bootstrap cụm
SEED <- 2026 # hạt giống ngẫu nhiên – giữ cố định để kết quả lặp lại được
# ---- Màu Okabe–Ito ----
MAU_NHOM <- c(Chung = "#999999", `Tân sinh` = "#D55E00", `Không tân sinh` = "#0072B2")
MAU_PP <- c(MBH = "#999999", BS = "#0072B2", CADx = "#E69F00")
# ---- Font chữ trong hình: Arial nếu máy có, nếu không dùng font không chân mặc định ----
FONT <- if ("Arial" %in% systemfonts::system_fonts()$family) "Arial" else "sans"
# ---- Định dạng số ----
fmtp <- function(p) ifelse(is.na(p), "—", ifelse(p < 0.001, "<0.001", sprintf("%.3f", p)))
nhan_p <- function(p) ifelse(is.na(p), "p = —", ifelse(p < 0.001, "p < 0.001", sprintf("p = %.3f", p))) # nhãn p trên hình
np <- function(x) { x <- as.numeric(x); pc <- if (sum(x) == 0) 0 * x else 100 * x / sum(x)
sprintf("%d (%.1f)", as.integer(x), pc) }
# ---- Định dạng bảng flextable thống nhất (Times New Roman 11, tiêu đề tô màu) ----
set_flextable_defaults(font.family = "Times New Roman", font.size = 11, padding = 3,
border.color = "#7F7F7F")
dinh_dang_bang <- function(ft, tieu_de = NULL, chu_thich = NULL) {
ft <- ft |> theme_box() |>
bg(part = "header", bg = "#DCE6F1") |> bold(part = "header") |>
align(part = "header", align = "center") |>
flextable::font(fontname = "Times New Roman", part = "all") |> fontsize(size = 11, part = "all") |>
autofit()
# Bảng rộng hơn khổ giấy (A4 dọc, lề trái 3 cm, phải 2 cm → 16 cm = 6.3 inch):
# thu hẹp cột theo tỷ lệ, GIỮ cỡ chữ 11, chữ tự xuống dòng trong ô
# (cột hẹp ≤ 0.8 inch như cột p giữ nguyên để "<0.001" không bị ngắt dòng)
rong <- dim(ft)$widths
if (sum(rong) > 6.3) {
hep <- rong <= 0.8
con_lai <- 6.3 - sum(rong[hep])
if (con_lai > 0 && any(!hep)) rong[!hep] <- rong[!hep] * con_lai / sum(rong[!hep]) else rong <- rong * 6.3 / sum(rong)
ft <- ft |> width(j = seq_along(rong), width = rong) |> set_table_properties(layout = "fixed")
}
ft <- ft
if (!is.null(chu_thich)) ft <- ft |> add_footer_lines(chu_thich) |> fontsize(size = 10, part = "footer") |>
italic(part = "footer")
if (!is.null(tieu_de)) ft <- ft |> set_caption(tieu_de)
ft
}
# ---- Giao diện biểu đồ chuẩn tạp chí ----
theme_bao <- function() {
ggpubr::theme_pubr(base_size = 9, base_family = FONT, legend = "top") +
theme(plot.title = element_blank(),
plot.subtitle = element_text(size = 8, colour = "grey25"),
plot.caption = element_text(size = 7, colour = "grey35", hjust = 0),
legend.title = element_blank(), legend.key.size = unit(3.5, "mm"),
axis.title = element_text(size = 9), axis.text = element_text(size = 8, colour = "black"),
plot.background = element_rect(fill = "white", colour = NA))
}
# ---- Lưu hình: PNG 300 dpi + TIFF 600 dpi nén LZW; kích thước theo cm ----
# rong_cm = 17.5 (2 cột tạp chí) hoặc 8.5 (1 cột)
luu_hinh <- function(p, ten_file, rong_cm = 17.5, cao_cm = 10) {
f_png <- file.path(thu_muc_kq, paste0(ten_file, ".png"))
f_tif <- file.path(thu_muc_kq, paste0(ten_file, ".tiff"))
ragg::agg_png(f_png, width = rong_cm, height = cao_cm, units = "cm", res = 300, background = "white")
print(p); invisible(dev.off())
ragg::agg_tiff(f_tif, width = rong_cm, height = cao_cm, units = "cm", res = 600,
background = "white", compression = "lzw")
print(p); invisible(dev.off())
invisible(c(f_png, f_tif))
}PHẦN 1. NHẬP DỮ LIỆU
1.1. Đọc file Excel
Đọc sheet “C” của file DATA_AI_POLYP.xlsx. Mỗi dòng là một polyp (đơn vị phân tích). Đọc tất cả cột dưới dạng văn bản (col_types = “text”) để giữ nguyên số 0 ở đầu mã bệnh nhân, sau đó chuyển Tuổi và điểm Boston về dạng số. Tên cột và nội dung được chuẩn hoá Unicode (NFC), vì Excel trên macOS đôi khi lưu chữ tiếng Việt có dấu theo kiểu tổ hợp khác.
raw <- read_excel(file_du_lieu, sheet = sheet_du_lieu, col_types = "text")
names(raw) <- stringi::stri_trans_nfc(trimws(names(raw)))
raw <- raw[rowSums(!is.na(raw)) > 0, ] # bỏ dòng trống hoàn toàn
cat("Số dòng (polyp):", nrow(raw), " | Số cột:", ncol(raw), "\n")## Số dòng (polyp): 275 | Số cột: 151.2. Xem nhanh dữ liệu
names(raw) # tên các cột## [1] "ID_POLYP" "Vị trí" "Paris"
## [4] "NICE" "Chẩn đoán bác sĩ" "Chẩn đoán CADx"
## [7] "Chẩn đoán mô bệnh học" "ID_BN" "Tuổi"
## [10] "Giới tính" "Lý do chỉ định" "Chuẩn bị ruột"
## [13] "Phân loại bác sĩ" "Phân loại CADx" "Phân loại mô bệnh học"head(raw, 5) # 5 dòng đầu## # A tibble: 5 × 15
## ID_POLYP `Vị trí` Paris NICE `Chẩn đoán bác sĩ` `Chẩn đoán CADx`
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 2001200460_POLYP 1 Đại tràng … Có c… Type… Adenoma Adenoma
## 2 2001200641_POLYP 1 Đại tràng … Lồi,… Type… Hyperplastic Hyperplastic
## 3 2000962527_POLYP 1 Đại tràng … Có c… Type… Adenoma Adenoma
## 4 0000101242_POLYP 1 Đại tràng … Lồi,… Type… Hyperplastic Hyperplastic
## 5 2000917107_POLYP 1 Đại tràng … Lồi,… Type… Hyperplastic Hyperplastic
## # ℹ 9 more variables: `Chẩn đoán mô bệnh học` <chr>, ID_BN <chr>, Tuổi <chr>,
## # `Giới tính` <chr>, `Lý do chỉ định` <chr>, `Chuẩn bị ruột` <chr>,
## # `Phân loại bác sĩ` <chr>, `Phân loại CADx` <chr>,
## # `Phân loại mô bệnh học` <chr>glimpse(raw) # kiểu dữ liệu từng cột## Rows: 275
## Columns: 15
## $ ID_POLYP <chr> "2001200460_POLYP 1", "2001200641_POLYP 1", "2…
## $ `Vị trí` <chr> "Đại tràng phải", "Đại tràng trái", "Đại tràng…
## $ Paris <chr> "Có cuống (0-Ip)", "Lồi, không cuống (0-Is)", …
## $ NICE <chr> "Type 2", "Type 1", "Type 2", "Type 1", "Type …
## $ `Chẩn đoán bác sĩ` <chr> "Adenoma", "Hyperplastic", "Adenoma", "Hyperpl…
## $ `Chẩn đoán CADx` <chr> "Adenoma", "Hyperplastic", "Adenoma", "Hyperpl…
## $ `Chẩn đoán mô bệnh học` <chr> "Adenoma", "Viêm", "Adenoma", "Hyperplastic", …
## $ ID_BN <chr> "2001200460", "2001200641", "2000962527", "000…
## $ Tuổi <chr> "71", "65", "73", "57", "26", "48", "39", "73"…
## $ `Giới tính` <chr> "Nữ", "Nữ", "Nam", "Nam", "Nữ", "Nam", "Nam", …
## $ `Lý do chỉ định` <chr> "Có triệu chứng", "Có triệu chứng", "Có triệu …
## $ `Chuẩn bị ruột` <chr> "8", "8", "9", "8", "7", "8", "7", "8", "8", "…
## $ `Phân loại bác sĩ` <chr> "Tân Sinh", "Không Tân Sinh", "Tân Sinh", "Khô…
## $ `Phân loại CADx` <chr> "Tân Sinh", "Không Tân Sinh", "Tân Sinh", "Khô…
## $ `Phân loại mô bệnh học` <chr> "Tân Sinh", "Không Tân Sinh", "Tân Sinh", "Khô…1.3. Đổi tên cột sang tên biến ngắn (không dấu)
ten_cot <- c("ID_POLYP" = "id_polyp", "Vị trí" = "vi_tri", "Paris" = "paris", "NICE" = "nice",
"Chẩn đoán bác sĩ" = "cd_bs", "Chẩn đoán CADx" = "cd_cadx",
"Chẩn đoán mô bệnh học" = "mbh", "ID_BN" = "id_bn", "Tuổi" = "tuoi",
"Giới tính" = "gioi", "Lý do chỉ định" = "ly_do", "Chuẩn bị ruột" = "boston",
"Phân loại bác sĩ" = "pl_bs", "Phân loại CADx" = "pl_cadx",
"Phân loại mô bệnh học" = "pl_mbh")
names(ten_cot) <- stringi::stri_trans_nfc(names(ten_cot))
thieu_cot <- setdiff(names(ten_cot), names(raw))
if (length(thieu_cot)) stop("File dữ liệu thiếu cột: ", paste(thieu_cot, collapse = ", "))
d0 <- raw[, names(ten_cot)]
names(d0) <- ten_cot
chuan_chu <- function(x) stringi::stri_trans_nfc(trimws(gsub("\\s+", " ", x)))
d0 <- d0 |> mutate(across(everything(), chuan_chu),
tuoi = as.numeric(tuoi), boston = as.numeric(boston))
str(d0)## tibble [275 × 15] (S3: tbl_df/tbl/data.frame)
## $ id_polyp: chr [1:275] "2001200460_POLYP 1" "2001200641_POLYP 1" "2000962527_POLYP 1" "0000101242_POLYP 1" ...
## $ vi_tri : chr [1:275] "Đại tràng phải" "Đại tràng trái" "Đại tràng Sigma" "Đại tràng trái" ...
## $ paris : chr [1:275] "Có cuống (0-Ip)" "Lồi, không cuống (0-Is)" "Có cuống (0-Ip)" "Lồi, không cuống (0-Is)" ...
## $ nice : chr [1:275] "Type 2" "Type 1" "Type 2" "Type 1" ...
## $ cd_bs : chr [1:275] "Adenoma" "Hyperplastic" "Adenoma" "Hyperplastic" ...
## $ cd_cadx : chr [1:275] "Adenoma" "Hyperplastic" "Adenoma" "Hyperplastic" ...
## $ mbh : chr [1:275] "Adenoma" "Viêm" "Adenoma" "Hyperplastic" ...
## $ id_bn : chr [1:275] "2001200460" "2001200641" "2000962527" "0000101242" ...
## $ tuoi : num [1:275] 71 65 73 57 26 48 39 73 54 40 ...
## $ gioi : chr [1:275] "Nữ" "Nữ" "Nam" "Nam" ...
## $ ly_do : chr [1:275] "Có triệu chứng" "Có triệu chứng" "Có triệu chứng" "Có triệu chứng" ...
## $ boston : num [1:275] 8 8 9 8 7 8 7 8 8 8 ...
## $ pl_bs : chr [1:275] "Tân Sinh" "Không Tân Sinh" "Tân Sinh" "Không Tân Sinh" ...
## $ pl_cadx : chr [1:275] "Tân Sinh" "Không Tân Sinh" "Tân Sinh" "Không Tân Sinh" ...
## $ pl_mbh : chr [1:275] "Tân Sinh" "Không Tân Sinh" "Tân Sinh" "Không Tân Sinh" ...#> tibble [275 × 15] (S3: tbl_df/tbl/data.frame)
#> $ id_polyp: chr [1:275] "2001200460_POLYP 1" "2001200641_POLYP 1" "2000962527_POLYP 1" "0000101242_POLYP 1" ...
#> $ vi_tri : chr [1:275] "Đại tràng phải" "Đại tràng trái" "Đại tràng Sigma" "Đại tràng trái" ...
#> $ paris : chr [1:275] "Có cuống (0-Ip)" "Lồi, không cuống (0-Is)" "Có cuống (0-Ip)" "Lồi, không cuống (0-Is)" ...
#> $ nice : chr [1:275] "Type 2" "Type 1" "Type 2" "Type 1" ...
#> $ cd_bs : chr [1:275] "Adenoma" "Hyperplastic" "Adenoma" "Hyperplastic" ...
#> $ cd_cadx : chr [1:275] "Adenoma" "Hyperplastic" "Adenoma" "Hyperplastic" ...
#> $ mbh : chr [1:275] "Adenoma" "Viêm" "Adenoma" "Hyperplastic" ...
#> $ id_bn : chr [1:275] "2001200460" "2001200641" "2000962527" "0000101242" ...
#> $ tuoi : num [1:275] 71 65 73 57 26 48 39 73 54 40 ...
#> $ gioi : chr [1:275] "Nữ" "Nữ" "Nam" "Nam" ...
#> $ ly_do : chr [1:275] "Có triệu chứng" "Có triệu chứng" "Có triệu chứng" "Có triệu chứng" ...
#> $ boston : num [1:275] 8 8 9 8 7 8 7 8 8 8 ...
#> $ pl_bs : chr [1:275] "Tân Sinh" "Không Tân Sinh" "Tân Sinh" "Không Tân Sinh" ...
#> $ pl_cadx : chr [1:275] "Tân Sinh" "Không Tân Sinh" "Tân Sinh" "Không Tân Sinh" ...
#> $ pl_mbh : chr [1:275] "Tân Sinh" "Không Tân Sinh" "Tân Sinh" "Không Tân Sinh" ...PHẦN 2. KIỂM TRA CHẤT LƯỢNG VÀ MÃ HOÁ BIẾN
2.1. Kiểm tra chất lượng dữ liệu (STARD mục 15–16)
tra <- list(
vi_tri = c("đại tràng phải" = 1, "đại tràng ngang" = 2, "đại tràng trái" = 3,
"đại tràng sigma" = 4, "trực tràng" = 5),
gioi = c("nam" = 1, "nữ" = 2),
ly_do = c("tầm soát" = 1, "có triệu chứng" = 2),
nice = c("type 1" = 1, "type 2" = 2, "type 3" = 3),
cd = c("adenoma" = 1, "ssl" = 2, "hyperplastic" = 3),
mbh = c("adenoma" = 1, "ssl" = 2, "hyperplastic" = 3, "carcinoma" = 4, "ung thư" = 4,
"viêm" = 5, "khác" = 5),
pl = c("tân sinh" = 1, "không tân sinh" = 0))
ma <- function(x, bang) unname(bang[tolower(x)])
# (a) Nhãn không hợp lệ
cap_kt <- list(c("vi_tri", "vi_tri"), c("gioi", "gioi"), c("ly_do", "ly_do"), c("nice", "nice"),
c("cd_bs", "cd"), c("cd_cadx", "cd"), c("mbh", "mbh"),
c("pl_bs", "pl"), c("pl_cadx", "pl"), c("pl_mbh", "pl"))
kt_nhan <- bind_rows(lapply(cap_kt, function(v) {
x <- d0[[v[1]]]; sai <- !is.na(x) & is.na(ma(x, tra[[v[2]]]))
if (any(sai)) data.frame(bien = v[1], id_polyp = d0$id_polyp[sai], gia_tri = x[sai]) else NULL }))
# (b) Mã polyp: trùng lặp, sai mẫu "mãBN_POLYP n", không khớp mã bệnh nhân
id_chuan <- toupper(gsub("[^A-Za-z0-9]", "", d0$id_polyp)) # bỏ khoảng trắng, gạch
trung <- d0$id_polyp[duplicated(id_chuan)]
sai_mau <- d0$id_polyp[!grepl("^[0-9]+_POLYP ?[0-9]+$", d0$id_polyp)]
lech_bn <- d0$id_polyp[!startsWith(id_chuan, toupper(d0$id_bn))]
# (c) Thông tin bệnh nhân không thống nhất giữa các polyp; giá trị ngoài khoảng
bn_mt <- d0 |> group_by(id_bn) |>
summarise(across(c(tuoi, gioi, ly_do, boston), n_distinct), .groups = "drop") |>
filter(tuoi > 1 | gioi > 1 | ly_do > 1 | boston > 1)
ngoai_khoang <- sum(d0$tuoi < 18 | d0$tuoi > 110 | d0$boston < 0 | d0$boston > 9, na.rm = TRUE)
# (d) Cột phân loại nhị phân có khớp với chẩn đoán chi tiết?
nhi_phan <- function(cd, bang) ifelse(is.na(cd), NA, as.integer(ma(cd, bang) %in% c(1, 2, 4)))
lech <- c(bs = sum(nhi_phan(d0$cd_bs, tra$cd) != ma(d0$pl_bs, tra$pl), na.rm = TRUE),
cadx = sum(nhi_phan(d0$cd_cadx, tra$cd) != ma(d0$pl_cadx, tra$pl), na.rm = TRUE),
mbh = sum(nhi_phan(d0$mbh, tra$mbh) != ma(d0$pl_mbh, tra$pl), na.rm = TRUE))
liet_ke <- function(v, toi_da = 5) if (!length(v)) "0" else
paste0(length(v), " (", paste(head(v, toi_da), collapse = "; "), if (length(v) > toi_da) "; …", ")")
kiem_tra <- data.frame(
`Nội dung kiểm tra` = c("Số polyp (dòng dữ liệu)", "Số bệnh nhân", "Nhãn không hợp lệ",
"Mã polyp trùng lặp (sau khi bỏ khoảng trắng, dấu gạch)",
"Mã polyp không đúng mẫu 'mãBN_POLYP n' hoặc 'mãBN_POLYPn' (chỉ cảnh báo)",
"Mã polyp không bắt đầu bằng mã bệnh nhân",
"Bệnh nhân có tuổi/giới/lý do/Boston khác nhau giữa các polyp",
"Tuổi hoặc điểm Boston ngoài khoảng hợp lệ",
"Thiếu chẩn đoán bác sĩ", "Thiếu chẩn đoán CADx", "Thiếu mô bệnh học",
"Phân loại nhị phân lệch với chẩn đoán chi tiết (BS / CADx / MBH)"),
`Kết quả` = c(nrow(d0), n_distinct(d0$id_bn), nrow(kt_nhan), liet_ke(unique(trung)),
liet_ke(sai_mau), liet_ke(lech_bn), nrow(bn_mt), ngoai_khoang,
sum(is.na(d0$cd_bs)), sum(is.na(d0$cd_cadx)), sum(is.na(d0$mbh)),
paste(lech, collapse = " / ")),
check.names = FALSE)
dinh_dang_bang(flextable(kiem_tra), "Bảng 0. Kiểm tra chất lượng dữ liệu")Nội dung kiểm tra | Kết quả |
|---|---|
Số polyp (dòng dữ liệu) | 275 |
Số bệnh nhân | 160 |
Nhãn không hợp lệ | 0 |
Mã polyp trùng lặp (sau khi bỏ khoảng trắng, dấu gạch) | 0 |
Mã polyp không đúng mẫu 'mãBN_POLYP n' hoặc 'mãBN_POLYPn' (chỉ cảnh báo) | 8 (2001203026- POLYP1; 2001072154_ POLYP1; 2001210640_ POLYP1; 2001045343_ POLYP1; 2001045957_POPLY1; …) |
Mã polyp không bắt đầu bằng mã bệnh nhân | 0 |
Bệnh nhân có tuổi/giới/lý do/Boston khác nhau giữa các polyp | 0 |
Tuổi hoặc điểm Boston ngoài khoảng hợp lệ | 0 |
Thiếu chẩn đoán bác sĩ | 0 |
Thiếu chẩn đoán CADx | 0 |
Thiếu mô bệnh học | 0 |
Phân loại nhị phân lệch với chẩn đoán chi tiết (BS / CADx / MBH) | 0 / 0 / 0 |
if (nrow(kt_nhan) > 0) {
print(kt_nhan)
stop("Có nhãn không hợp lệ – sửa file dữ liệu theo danh sách trên rồi chạy lại.")
}2.2. Mã hoá biến phân tích
Tạo 3 bộ dữ liệu: d (tất cả polyp), d_pair (polyp có đủ chẩn đoán bác sĩ, CADx và mô bệnh học – dùng cho phân tích bắt cặp), bn (mỗi dòng một bệnh nhân – chỉ dùng cho Bảng 3.1).
lab_vt <- c("Đại tràng phải", "Đại tràng ngang", "Đại tràng trái", "Đại tràng sigma", "Trực tràng")
lab_mbh <- c("U tuyến", "SSL", "Tăng sản", "Carcinoma", "Viêm/khác")
lv3 <- c("U tuyến", "SSL", "Tăng sản/không tân sinh")
paris_lv <- c("Có cuống (0-Ip)", "Lồi, không cuống (0-Is)", "Bán cuống (0-Isp)",
"Phẳng gồ (0-IIa)", "Phẳng (0-IIb)", "Phẳng lõm (0-IIc)")
paris_la <- setdiff(unique(d0$paris), paris_lv)
d <- d0 |>
mutate(
vi_tri_n = ma(vi_tri, tra$vi_tri),
vi_tri_f = factor(lab_vt[vi_tri_n], levels = lab_vt),
paris_f = droplevels(factor(paris, levels = c(paris_lv, paris_la))),
nice_f = droplevels(factor(paste("NICE", tolower(nice)), levels = paste("NICE type", 1:3))),
gioi_f = factor(c("Nam", "Nữ")[ma(gioi, tra$gioi)], levels = c("Nam", "Nữ")),
ly_do_f = factor(c("Tầm soát", "Có triệu chứng")[ma(ly_do, tra$ly_do)],
levels = c("Tầm soát", "Có triệu chứng")),
bs_n = ma(cd_bs, tra$cd), cadx_n = ma(cd_cadx, tra$cd), mbh_n = ma(mbh, tra$mbh),
mbh_f = factor(lab_mbh[mbh_n], levels = lab_mbh),
ref_ts = as.integer(mbh_n %in% c(1, 2, 4)),
bs_ts = ifelse(is.na(bs_n), NA_integer_, as.integer(bs_n %in% c(1, 2))),
cadx_ts = ifelse(is.na(cadx_n), NA_integer_, as.integer(cadx_n %in% c(1, 2))),
nhom = factor(ifelse(ref_ts == 1, "Tân sinh", "Không tân sinh"),
levels = c("Tân sinh", "Không tân sinh")),
ref3 = factor(lv3[c(1, 2, 3, 1, 3)[mbh_n]], levels = lv3),
bs3 = factor(lv3[bs_n], levels = lv3),
cadx3 = factor(lv3[cadx_n], levels = lv3),
vung = factor(ifelse(vi_tri_n %in% c(4, 5), "Trực tràng – sigma", "Đại tràng phải, ngang, trái"),
levels = c("Trực tràng – sigma", "Đại tràng phải, ngang, trái")))
d_pair <- d |> filter(!is.na(mbh_n), !is.na(bs_ts), !is.na(cadx_ts))
bn <- d |> group_by(id_bn) |>
summarise(tuoi = first(tuoi), gioi_f = first(gioi_f), ly_do_f = first(ly_do_f),
boston = first(boston), so_polyp = n(), co_ts = any(ref_ts == 1), .groups = "drop") |>
mutate(nhom = factor(ifelse(co_ts, "Tân sinh", "Không tân sinh"),
levels = c("Tân sinh", "Không tân sinh")))
cat(sprintf("Polyp: %d | Polyp phân tích bắt cặp: %d (tân sinh %d, không tân sinh %d) | Bệnh nhân: %d\n",
nrow(d), nrow(d_pair), sum(d_pair$ref_ts == 1), sum(d_pair$ref_ts == 0), nrow(bn)))## Polyp: 275 | Polyp phân tích bắt cặp: 275 (tân sinh 130, không tân sinh 145) | Bệnh nhân: 1602.3. Các hàm thống kê dùng chung
Các hàm dưới đây được dùng lại ở nhiều bảng:
chi_so(): độ nhạy (Se), độ đặc hiệu (Sp), PPV, NPV, độ chính xác (Acc), LR+, LR−.
boot_cum(): bootstrap cụm – lấy mẫu lại theo bệnh nhân (không theo polyp), vì các polyp của cùng một bệnh nhân không độc lập. Từ đó tính KTC 95% bách phân vị và p bootstrap.
mcnemar_cx(): kiểm định McNemar chính xác (phân phối nhị thức trên các cặp bất đồng).
gee_p(): GEE logit, tương quan trao đổi (exchangeable), cụm = bệnh nhân – so sánh tỷ lệ “chẩn đoán đúng” giữa bác sĩ và CADx, có hiệu chỉnh cụm.
kappa_z(): hệ số Kappa Cohen và kiểm định Kappa = 0 (phương sai dưới H0 theo Fleiss, Cohen & Everitt 1969).
bowker(): kiểm định Bowker (đối xứng) cho bảng bắt cặp 3 × 3.
p_phan_loai(), p_lien_tuc(): kiểm định so sánh 2 nhóm, viết theo đúng quy tắc của compareGroups để bảng table1 và compareGroups cho cùng p-value.
# ---- Chỉ số chẩn đoán ----
dem <- function(test, ref) c(TP = sum(test == 1 & ref == 1), FP = sum(test == 1 & ref == 0),
FN = sum(test == 0 & ref == 1), TN = sum(test == 0 & ref == 0))
chi_so <- function(test, ref) {
k <- dem(test, ref); TP <- k[["TP"]]; FP <- k[["FP"]]; FN <- k[["FN"]]; TN <- k[["TN"]]
se <- TP / (TP + FN); sp <- TN / (TN + FP)
c(Se = se, Sp = sp, PPV = TP / (TP + FP), NPV = TN / (TN + FN), Acc = (TP + TN) / sum(k),
LRp = se / (1 - sp), LRn = (1 - se) / sp)
}
# ---- Bootstrap cụm theo bệnh nhân ----
boot_cum <- function(data, stat_fun, B, seed) {
set.seed(seed); ds <- split(data, data$id_bn); ids <- names(ds)
t(replicate(B, stat_fun(bind_rows(ds[sample(ids, length(ids), replace = TRUE)]))))
}
ci_q <- function(m) apply(m, 2, function(v) quantile(v[is.finite(v)], c(.025, .975), na.rm = TRUE))
p_boot <- function(v) { v <- v[is.finite(v)]; min(1, 2 * min(mean(v <= 0), mean(v >= 0))) }
# ---- KTC Wilson (không hiệu chỉnh cụm) ----
wilson <- function(x, n, z = qnorm(.975)) {
if (is.na(n) || n == 0) return(c(NA, NA, NA)); p <- x / n
c0 <- (p + z^2 / (2 * n)) / (1 + z^2 / n); h <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2)) / (1 + z^2 / n)
c(p, max(0, c0 - h), min(1, c0 + h))
}
fmt_ci <- function(v, k = 1) ifelse(is.na(v[1]), "—",
sprintf("%.*f (%.*f–%.*f)", k, 100 * v[1], k, 100 * v[2], k, 100 * v[3]))
# ---- McNemar chính xác: a, b = véc-tơ TRUE/FALSE (đúng/sai hoặc có/không) trên cùng polyp ----
mcnemar_cx <- function(a, b) { n10 <- sum(a & !b); n01 <- sum(!a & b)
if (n10 + n01 == 0) NA else binom.test(min(n10, n01), n10 + n01)$p.value }
# ---- GEE: so sánh tỷ lệ chẩn đoán đúng BS vs CADx, cụm = bệnh nhân ----
gee_p <- function(sub) {
long <- bind_rows(transmute(sub, id_bn, pp = "BS", dung = as.integer(bs_ts == ref_ts)),
transmute(sub, id_bn, pp = "CADx", dung = as.integer(cadx_ts == ref_ts))) |>
arrange(id_bn) |> mutate(pp = factor(pp, levels = c("BS", "CADx")), id = as.integer(factor(id_bn)))
m <- try(geeglm(dung ~ pp, family = binomial, id = id, data = long, corstr = "exchangeable"), silent = TRUE)
if (inherits(m, "try-error")) NA else summary(m)$coefficients["ppCADx", "Pr(>|W|)"]
}
# ---- Kappa Cohen + kiểm định Kappa = 0 ----
kappa_z <- function(a, b, lv) {
t <- table(factor(a, lv), factor(b, lv)); n <- sum(t); p <- t / n
pr <- rowSums(p); pc <- colSums(p); po <- sum(diag(p)); pe <- sum(pr * pc)
k <- (po - pe) / (1 - pe)
se0 <- sqrt(pe + pe^2 - sum(pr * pc * (pr + pc))) / ((1 - pe) * sqrt(n))
c(po = po, k = k, z = k / se0, p = 2 * pnorm(-abs(k / se0)))
}
# ---- Bowker (đối xứng) cho bảng vuông k × k ----
bowker <- function(a, b, lv) {
t <- table(factor(a, lv), factor(b, lv)); st <- 0; df <- 0
for (i in 1:(length(lv) - 1)) for (j in (i + 1):length(lv)) {
s <- t[i, j] + t[j, i]; if (s > 0) { st <- st + (t[i, j] - t[j, i])^2 / s; df <- df + 1 } }
if (df == 0) return(c(chi2 = NA, df = 0, p = NA))
c(chi2 = st, df = df, p = pchisq(st, df, lower.tail = FALSE))
}
# ---- Kiểm định 2 nhóm theo đúng quy tắc compareGroups ----
p_phan_loai <- function(x, g) { # χ² (mặc định R); Fisher chính xác nếu tần số mong đợi < 5
t <- table(x, g); t <- t[rowSums(t) > 0, colSums(t) > 0, drop = FALSE]
if (nrow(t) < 2 || ncol(t) < 2) return(NA)
e <- outer(rowSums(t), colSums(t)) / sum(t)
if (any(e < 5)) fisher.test(t)$p.value else chisq.test(t)$p.value
}
p_lien_tuc <- function(x, g, phan_phoi_chuan) { # t-test Welch nếu chuẩn; Mann-Whitney nếu không
if (phan_phoi_chuan) t.test(x ~ g, var.equal = FALSE)$p.value else kruskal.test(x ~ g)$p.value
}PHẦN 3. PHÂN TÍCH ĐẶC ĐIỂM MẪU NGHIÊN CỨU (STARD mục 20)
** 3.1. Bảng 3.1 – Đặc điểm bệnh nhân (đơn vị: bệnh nhân)**
Đơn vị phân tích là bệnh nhân (n = 160), để mỗi người chỉ được đếm một lần và p-value hợp lệ.
Chọn kiểm định (compareGroups tự động):
Biến liên tục (Tuổi, Số polyp/bệnh nhân): kiểm tra phân phối chuẩn bằng Shapiro–Wilk (α = 0,05). Nếu chuẩn → trình bày TB ± ĐLC, dùng t-test Welch; nếu không → trình bày trung vị [Q1; Q3], dùng Mann–Whitney.
Điểm Boston (7–9) là biến thứ bậc: trình bày cả dạng số (trung vị [Q1; Q3], Mann–Whitney) và dạng phân loại (n (%), χ²/Fisher).
Biến phân loại: χ²; dùng Fisher chính xác khi có ô tần số mong đợi < 5.
Bước (a) gắn nhãn tiếng Việt cho biến, (b) chạy compareGroups(), (c) tạo bảng bằng createTable(), (d) chuyển thành bảng Word, (e) xuất riêng file Word/Excel của compareGroups.
# (a) Dữ liệu mức bệnh nhân + nhãn biến tiếng Việt
x31 <- bn |>
mutate(nhom_tuoi40 = factor(ifelse(tuoi < 40, "< 40", "≥ 40"), levels = c("< 40", "≥ 40")),
nhom_tuoi50 = factor(ifelse(tuoi < 50, "< 50", "≥ 50"), levels = c("< 50", "≥ 50")),
boston_f = factor(boston, levels = sort(unique(boston))))
nhan31 <- c(tuoi = "Tuổi (năm)", nhom_tuoi40 = "Nhóm tuổi (mốc 40)", nhom_tuoi50 = "Nhóm tuổi (mốc 50)",
gioi_f = "Giới", ly_do_f = "Lý do nội soi", boston = "Điểm Boston (dạng số)",
boston_f = "Điểm Boston (phân loại)", so_polyp = "Số polyp ≤ 5 mm/bệnh nhân")
for (v in names(nhan31)) attr(x31[[v]], "label") <- nhan31[[v]]
# Kiểm tra phân phối chuẩn (Shapiro–Wilk) – quyết định cách trình bày biến liên tục
chuan31 <- c(tuoi = shapiro.test(x31$tuoi)$p.value > 0.05,
so_polyp = shapiro.test(x31$so_polyp)$p.value > 0.05,
boston = FALSE) # Boston: biến thứ bậc → luôn phi tham số
chuan31## tuoi so_polyp boston
## TRUE FALSE FALSE# (b) compareGroups: method = 1 (chuẩn), 2 (không chuẩn), NA (tự kiểm tra Shapiro–Wilk)
cg31 <- compareGroups(nhom ~ tuoi + nhom_tuoi40 + nhom_tuoi50 + gioi_f + ly_do_f +
boston + boston_f + so_polyp,
data = x31, method = c(tuoi = NA, so_polyp = NA, boston = 2),
var.equal = FALSE) # t-test Welch
# (c) createTable: show.all = cột Chung; sd.type = 2 → "TB ± ĐLC"; digits = 1 chữ số thập phân
tab31 <- createTable(cg31, show.all = TRUE, sd.type = 2, digits = 1, show.p.overall = TRUE)
print(tab31) # xem nhanh trong Console##
## --------Summary descriptives table by 'nhom'---------
##
## ______________________________________________________________________________
## [ALL] Tân sinh Không tân sinh p.overall
## N=160 N=90 N=70
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Tuổi (năm) 54.6±12.8 56.0±13.4 52.9±11.9 0.135
## Nhóm tuổi (mốc 40): 1.000
## < 40 26 (16.2%) 15 (16.7%) 11 (15.7%)
## ≥ 40 134 (83.8%) 75 (83.3%) 59 (84.3%)
## Nhóm tuổi (mốc 50): 0.949
## < 50 51 (31.9%) 28 (31.1%) 23 (32.9%)
## ≥ 50 109 (68.1%) 62 (68.9%) 47 (67.1%)
## Giới: 0.626
## Nam 64 (40.0%) 34 (37.8%) 30 (42.9%)
## Nữ 96 (60.0%) 56 (62.2%) 40 (57.1%)
## Lý do nội soi: 0.608
## Tầm soát 15 (9.4%) 7 (7.8%) 8 (11.4%)
## Có triệu chứng 145 (90.6%) 83 (92.2%) 62 (88.6%)
## Điểm Boston (dạng số) 8.0 [8.0;9.0] 8.0 [8.0;9.0] 8.0 [8.0;8.0] 0.032
## Điểm Boston (phân loại): 0.096
## 7 20 (12.5%) 8 (8.9%) 12 (17.1%)
## 8 99 (61.9%) 54 (60.0%) 45 (64.3%)
## 9 41 (25.6%) 28 (31.1%) 13 (18.6%)
## Số polyp ≤ 5 mm/bệnh nhân 1.0 [1.0;2.0] 2.0 [1.0;3.0] 1.0 [1.0;2.0] <0.001
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯# Hàm chuyển kết quả createTable() thành data.frame 5 cột: Đặc điểm | Chung | Tân sinh | Không tân sinh | p
cg_sang_bang <- function(tab) {
m <- tab$descr; rn <- rownames(m); out <- list(); da_co <- character(0)
for (i in seq_along(rn)) {
if (grepl(": ", rn[i], fixed = TRUE)) { # dòng thuộc biến phân loại "Biến: mức"
bien <- sub(": .*$", "", rn[i]); muc <- sub("^[^:]*: ", "", rn[i])
if (!bien %in% da_co) { # dòng tiêu đề biến, mang p-value
out[[length(out) + 1]] <- c(bien, "", "", "", m[i, 4]); da_co <- c(da_co, bien) }
out[[length(out) + 1]] <- c(paste0(" ", muc), m[i, 1], m[i, 2], m[i, 3], "")
} else out[[length(out) + 1]] <- c(rn[i], m[i, 1], m[i, 2], m[i, 3], m[i, 4])
}
df <- as.data.frame(do.call(rbind, out), stringsAsFactors = FALSE)
names(df) <- c("Đặc điểm", "Chung", "Tân sinh", "Không tân sinh", "p")
df$p[is.na(df$p)] <- ""
df
}
# (d) Bảng Word# (d) Bảng Word
tb31 <- cg_sang_bang(tab31)
n31 <- table(x31$nhom)
ft31 <- flextable(tb31) |>
set_header_labels(Chung = sprintf("Chung\n(n = %d)", nrow(x31)),
`Tân sinh` = sprintf("Tân sinh\n(n = %d)", n31[["Tân sinh"]]),
`Không tân sinh` = sprintf("Không tân sinh\n(n = %d)", n31[["Không tân sinh"]])) |>
align(j = 2:5, align = "center", part = "body")
ft31 <- dinh_dang_bang(ft31, "Bảng 3.1. Đặc điểm của bệnh nhân trong nghiên cứu theo nhóm mô bệnh học",
c("Đơn vị: bệnh nhân. Nhóm Tân sinh: bệnh nhân có ≥ 1 polyp tân sinh trên mô bệnh học; nhóm Không tân sinh: chỉ có polyp không tân sinh.",
"Số liệu: n (%), TB ± ĐLC hoặc trung vị [Q1; Q3] (chọn theo kiểm định Shapiro–Wilk).",
"p: t-test Welch (TB ± ĐLC); Mann–Whitney (trung vị); χ² hoặc Fisher chính xác (tỷ lệ; Fisher khi tần số mong đợi < 5)."))
ft31Đặc điểm | Chung | Tân sinh | Không tân sinh | p |
|---|---|---|---|---|
Tuổi (năm) | 54.6±12.8 | 56.0±13.4 | 52.9±11.9 | 0.135 |
Nhóm tuổi (mốc 40) | 1.000 | |||
< 40 | 26 (16.2%) | 15 (16.7%) | 11 (15.7%) | |
≥ 40 | 134 (83.8%) | 75 (83.3%) | 59 (84.3%) | |
Nhóm tuổi (mốc 50) | 0.949 | |||
< 50 | 51 (31.9%) | 28 (31.1%) | 23 (32.9%) | |
≥ 50 | 109 (68.1%) | 62 (68.9%) | 47 (67.1%) | |
Giới | 0.626 | |||
Nam | 64 (40.0%) | 34 (37.8%) | 30 (42.9%) | |
Nữ | 96 (60.0%) | 56 (62.2%) | 40 (57.1%) | |
Lý do nội soi | 0.608 | |||
Tầm soát | 15 (9.4%) | 7 (7.8%) | 8 (11.4%) | |
Có triệu chứng | 145 (90.6%) | 83 (92.2%) | 62 (88.6%) | |
Điểm Boston (dạng số) | 8.0 [8.0;9.0] | 8.0 [8.0;9.0] | 8.0 [8.0;8.0] | 0.032 |
Điểm Boston (phân loại) | 0.096 | |||
7 | 20 (12.5%) | 8 (8.9%) | 12 (17.1%) | |
8 | 99 (61.9%) | 54 (60.0%) | 45 (64.3%) | |
9 | 41 (25.6%) | 28 (31.1%) | 13 (18.6%) | |
Số polyp ≤ 5 mm/bệnh nhân | 1.0 [1.0;2.0] | 2.0 [1.0;3.0] | 1.0 [1.0;2.0] | <0.001 |
Đơn vị: bệnh nhân. Nhóm Tân sinh: bệnh nhân có ≥ 1 polyp tân sinh trên mô bệnh học; nhóm Không tân sinh: chỉ có polyp không tân sinh. | ||||
Số liệu: n (%), TB ± ĐLC hoặc trung vị [Q1; Q3] (chọn theo kiểm định Shapiro–Wilk). | ||||
p: t-test Welch (TB ± ĐLC); Mann–Whitney (trung vị); χ² hoặc Fisher chính xác (tỷ lệ; Fisher khi tần số mong đợi < 5). | ||||
# (e) Xuất riêng Bảng 3.1: Word (flextable, TNR 11) và Excel (compareGroups)
# Ghi chú: không dùng export2word() của compareGroups vì hàm này tự gọi rmarkdown bên trong,
# gây lỗi khi đang Knit; save_as_docx() cho kết quả giống bảng trong báo cáo.
save_as_docx(ft31, path = file.path(thu_muc_kq, "Bang_3_1_compareGroups.docx"))
export2xls(tab31, file = file.path(thu_muc_kq, "Bang_3_1_compareGroups.xlsx"))Bảng 3.1 bằng table1 (bản HTML, có cột Chung và p-value)
table1 không tự tính p-value, nên viết hàm p_table1() gắn vào tham số extra.col. Hàm này dùng đúng quy tắc của compareGroups (Shapiro–Wilk → t-test Welch/Mann–Whitney; χ²/Fisher), vì vậy p-value ở hai bảng phải trùng nhau – khối lệnh cuối mục này đối chiếu tự động.
# Hàm hiển thị: biến liên tục theo kết quả Shapiro–Wilk; biến phân loại "n (%)"
hien_thi_t1 <- function(x, name, ...) {
if (is.numeric(x)) {
if (isTRUE(chuan31[name])) c("", "TB ± ĐLC" = sprintf("%.1f ± %.1f", mean(x), sd(x)))
else { q <- quantile(x, c(.25, .5, .75))
c("", "Trung vị [Q1; Q3]" = sprintf("%.1f [%.1f; %.1f]", q[2], q[1], q[3])) }
} else { t <- table(x); c("", setNames(sprintf("%d (%.1f%%)", t, 100 * t / sum(t)), names(t))) }
}
# Hàm p-value cho table1 (x = danh sách giá trị theo từng nhóm)
p_table1 <- function(x, name, ...) {
x <- x[names(x) != "Chung"] # bỏ cột Chung (overall), chỉ so sánh 2 nhóm
y <- unlist(x); g <- factor(rep(seq_along(x), times = sapply(x, length)))
p <- if (is.numeric(y)) p_lien_tuc(y, g, isTRUE(chuan31[name])) else p_phan_loai(y, g)
c("", fmtp(p))
}
bang31_t1 <- table1(~ tuoi + nhom_tuoi40 + nhom_tuoi50 + gioi_f + ly_do_f + boston + boston_f + so_polyp | nhom,
data = x31, overall = "Chung", render = hien_thi_t1,
extra.col = list(p = p_table1), extra.col.pos = 4,
caption = "Bảng 3.1 (table1). Đặc điểm bệnh nhân theo nhóm mô bệnh học")
bang31_t1| Tân sinh (N=90) |
Không tân sinh (N=70) |
Chung (N=160) |
p | |
|---|---|---|---|---|
| Tuổi (năm) | ||||
| TB ± ĐLC | 56.0 ± 13.4 | 52.9 ± 11.9 | 54.6 ± 12.8 | 0.135 |
| Nhóm tuổi (mốc 40) | ||||
| < 40 | 15 (16.7%) | 11 (15.7%) | 26 (16.2%) | 1.000 |
| ≥ 40 | 75 (83.3%) | 59 (84.3%) | 134 (83.8%) | |
| Nhóm tuổi (mốc 50) | ||||
| < 50 | 28 (31.1%) | 23 (32.9%) | 51 (31.9%) | 0.949 |
| ≥ 50 | 62 (68.9%) | 47 (67.1%) | 109 (68.1%) | |
| Giới | ||||
| Nam | 34 (37.8%) | 30 (42.9%) | 64 (40.0%) | 0.626 |
| Nữ | 56 (62.2%) | 40 (57.1%) | 96 (60.0%) | |
| Lý do nội soi | ||||
| Tầm soát | 7 (7.8%) | 8 (11.4%) | 15 (9.4%) | 0.608 |
| Có triệu chứng | 83 (92.2%) | 62 (88.6%) | 145 (90.6%) | |
| Điểm Boston (dạng số) | ||||
| Trung vị [Q1; Q3] | 8.0 [8.0; 9.0] | 8.0 [8.0; 8.0] | 8.0 [8.0; 9.0] | 0.032 |
| Điểm Boston (phân loại) | ||||
| 7 | 8 (8.9%) | 12 (17.1%) | 20 (12.5%) | 0.096 |
| 8 | 54 (60.0%) | 45 (64.3%) | 99 (61.9%) | |
| 9 | 28 (31.1%) | 13 (18.6%) | 41 (25.6%) | |
| Số polyp ≤ 5 mm/bệnh nhân | ||||
| Trung vị [Q1; Q3] | 2.0 [1.0; 3.0] | 1.0 [1.0; 2.0] | 1.0 [1.0; 2.0] | <0.001 |
# Đối chiếu p-value giữa compareGroups và cách tính trong table1
p_t1 <- sapply(names(nhan31), function(v) {
x <- x31[[v]]; if (is.numeric(x)) p_lien_tuc(x, x31$nhom, isTRUE(chuan31[v])) else p_phan_loai(x, x31$nhom) })
p_cg <- tb31$p[tb31$p != ""]
doi_chieu <- data.frame(`Biến` = nhan31, `p compareGroups` = p_cg, `p table1` = fmtp(p_t1),
`Khớp` = ifelse(p_cg == fmtp(p_t1), "✔", "✘ kiểm tra"), check.names = FALSE)
knitr::kable(doi_chieu, row.names = FALSE, caption = "Đối chiếu p-value giữa compareGroups và table1")| Biến | p compareGroups | p table1 | Khớp |
|---|---|---|---|
| Tuổi (năm) | 0.135 | 0.135 | ✔ |
| Nhóm tuổi (mốc 40) | 1.000 | 1.000 | ✔ |
| Nhóm tuổi (mốc 50) | 0.949 | 0.949 | ✔ |
| Giới | 0.626 | 0.626 | ✔ |
| Lý do nội soi | 0.608 | 0.608 | ✔ |
| Điểm Boston (dạng số) | 0.032 | 0.032 | ✔ |
| Điểm Boston (phân loại) | 0.096 | 0.096 | ✔ |
| Số polyp ≤ 5 mm/bệnh nhân | <0.001 | <0.001 | ✔ |
3.2. Bảng 3.2 – Đặc điểm polyp (đơn vị: polyp)
Từ bảng này trở đi đơn vị phân tích là polyp (n = 275). Vì một bệnh nhân có thể có nhiều polyp, ngoài p của χ²/Fisher (compareGroups) còn thêm cột p hiệu chỉnh cụm (GEE): mô hình logit tân sinh ~ đặc điểm, cụm = bệnh nhân, tương quan trao đổi, kiểm định Wald.
x32 <- d |> select(id_bn, ref_ts, nhom, vi_tri_f, paris_f, nice_f)
nhan32 <- c(vi_tri_f = "Vị trí", paris_f = "Hình thái (Paris)", nice_f = "Phân loại NICE")
for (v in names(nhan32)) attr(x32[[v]], "label") <- nhan32[[v]]
cg32 <- compareGroups(nhom ~ vi_tri_f + paris_f + nice_f, data = x32)
tab32 <- createTable(cg32, show.all = TRUE, digits = 1, show.p.overall = TRUE)
export2xls(tab32, file = file.path(thu_muc_kq, "Bang_3_2_compareGroups.xlsx"))
# p hiệu chỉnh cụm theo bệnh nhân (GEE)
gee_bien <- function(bien) {
dd <- x32 |> arrange(id_bn) |> mutate(id = as.integer(factor(id_bn)), x = .data[[bien]])
m <- geeglm(ref_ts ~ x, family = binomial, id = id, data = dd, corstr = "exchangeable")
anova(m)[["P(>|Chi|)"]][1]
}
p_gee32 <- sapply(names(nhan32), gee_bien)
tb32 <- cg_sang_bang(tab32)
tb32$`p (GEE)` <- ""
tb32$`p (GEE)`[match(nhan32, tb32$`Đặc điểm`)] <- fmtp(p_gee32)
names(tb32)[5] <- "p (χ²/Fisher)"
ft32 <- flextable(tb32) |>
set_header_labels(Chung = sprintf("Chung\n(n = %d)", nrow(d)),
`Tân sinh` = sprintf("Tân sinh\n(n = %d)", sum(d$nhom == "Tân sinh")),
`Không tân sinh` = sprintf("Không tân sinh\n(n = %d)", sum(d$nhom == "Không tân sinh"))) |>
align(j = 2:6, align = "center", part = "body")
ft32 <- dinh_dang_bang(ft32, "Bảng 3.2. Đặc điểm polyp theo nhóm mô bệnh học",
c("Đơn vị: polyp; nhóm theo mô bệnh học của từng polyp. Số liệu: n (%).",
"p (χ²/Fisher): không hiệu chỉnh cụm; p (GEE): hồi quy logistic GEE, cụm = bệnh nhân, tương quan trao đổi, kiểm định Wald."))
ft32Đặc điểm | Chung | Tân sinh | Không tân sinh | p (χ²/Fisher) | p (GEE) |
|---|---|---|---|---|---|
Vị trí | <0.001 | <0.001 | |||
Đại tràng phải | 61 (22.2%) | 36 (27.7%) | 25 (17.2%) | ||
Đại tràng ngang | 48 (17.5%) | 31 (23.8%) | 17 (11.7%) | ||
Đại tràng trái | 54 (19.6%) | 28 (21.5%) | 26 (17.9%) | ||
Đại tràng sigma | 66 (24.0%) | 27 (20.8%) | 39 (26.9%) | ||
Trực tràng | 46 (16.7%) | 8 (6.2%) | 38 (26.2%) | ||
Hình thái (Paris) | <0.001 | <0.001 | |||
Có cuống (0-Ip) | 191 (69.5%) | 111 (85.4%) | 80 (55.2%) | ||
Lồi, không cuống (0-Is) | 84 (30.5%) | 19 (14.6%) | 65 (44.8%) | ||
Phân loại NICE | <0.001 | <0.001 | |||
NICE type 1 | 145 (52.7%) | 39 (30.0%) | 106 (73.1%) | ||
NICE type 2 | 130 (47.3%) | 91 (70.0%) | 39 (26.9%) | ||
Đơn vị: polyp; nhóm theo mô bệnh học của từng polyp. Số liệu: n (%). | |||||
p (χ²/Fisher): không hiệu chỉnh cụm; p (GEE): hồi quy logistic GEE, cụm = bệnh nhân, tương quan trao đổi, kiểm định Wald. | |||||
save_as_docx(ft32, path = file.path(thu_muc_kq, "Bang_3_2_compareGroups.docx"))table1(~ vi_tri_f + paris_f + nice_f | nhom, data = x32, overall = "Chung", render = hien_thi_t1,
extra.col = list(`p (χ²/Fisher)` = p_table1), extra.col.pos = 4,
caption = "Bảng 3.2 (table1). Đặc điểm polyp theo nhóm mô bệnh học")| Tân sinh (N=130) |
Không tân sinh (N=145) |
Chung (N=275) |
p (χ²/Fisher) | |
|---|---|---|---|---|
| Vị trí | ||||
| Đại tràng phải | 36 (27.7%) | 25 (17.2%) | 61 (22.2%) | <0.001 |
| Đại tràng ngang | 31 (23.8%) | 17 (11.7%) | 48 (17.5%) | |
| Đại tràng trái | 28 (21.5%) | 26 (17.9%) | 54 (19.6%) | |
| Đại tràng sigma | 27 (20.8%) | 39 (26.9%) | 66 (24.0%) | |
| Trực tràng | 8 (6.2%) | 38 (26.2%) | 46 (16.7%) | |
| Hình thái (Paris) | ||||
| Có cuống (0-Ip) | 111 (85.4%) | 80 (55.2%) | 191 (69.5%) | <0.001 |
| Lồi, không cuống (0-Is) | 19 (14.6%) | 65 (44.8%) | 84 (30.5%) | |
| Phân loại NICE | ||||
| NICE type 1 | 39 (30.0%) | 106 (73.1%) | 145 (52.7%) | <0.001 |
| NICE type 2 | 91 (70.0%) | 39 (26.9%) | 130 (47.3%) |
3.3. Biểu đồ 3.1 – Vị trí polyp theo nhóm mô bệnh học
Mỗi biểu đồ được viết thành một hàm có tham số ngôn ngữ (“vi” hoặc “en”), gọi hai lần để có bản tiếng Việt (luận văn) và bản tiếng Anh (nộp tạp chí). Mỗi bản được lưu thành PNG 300 dpi và TIFF 600 dpi (LZW), khổ 2 cột tạp chí (17,5 cm). p-value ghi ở phụ đề.
p31_chi <- p_phan_loai(d$vi_tri_f, d$nhom); p31_gee <- p_gee32[["vi_tri_f"]]
vt <- bind_rows(d |> count(vi_tri_f) |> mutate(nhom = "Chung"),
d |> count(nhom, vi_tri_f) |> mutate(nhom = as.character(nhom))) |>
group_by(nhom) |> mutate(tl = n / sum(n)) |> ungroup() |>
mutate(nhom = factor(nhom, levels = names(MAU_NHOM)))
ve_31 <- function(lang = "vi") {
en <- lang == "en"
nhan_vt <- if (en) c("Right colon", "Transverse colon", "Left colon", "Sigmoid colon", "Rectum") else lab_vt
nhan_nh <- if (en) c(sprintf("Overall (n = %d)", nrow(d)), sprintf("Neoplastic (n = %d)", sum(d$nhom == "Tân sinh")),
sprintf("Non-neoplastic (n = %d)", sum(d$nhom == "Không tân sinh")))
else c(sprintf("Chung (n = %d)", nrow(d)), sprintf("Tân sinh (n = %d)", sum(d$nhom == "Tân sinh")),
sprintf("Không tân sinh (n = %d)", sum(d$nhom == "Không tân sinh")))
ggplot(vt, aes(vi_tri_f, tl, fill = nhom)) +
geom_col(position = position_dodge(0.85), width = 0.8, colour = "white", linewidth = 0.2) +
geom_text(aes(label = sprintf("%d\n(%.1f%%)", n, 100 * tl)), position = position_dodge(0.85),
vjust = -0.15, size = 2.3, lineheight = 0.85, family = FONT) +
scale_fill_manual(values = MAU_NHOM, labels = nhan_nh) +
scale_x_discrete(labels = nhan_vt) +
scale_y_continuous(labels = percent, expand = expansion(mult = c(0, 0.18))) +
labs(x = NULL, y = if (en) "Proportion of polyps within group" else "Tỷ lệ polyp trong nhóm",
subtitle = if (en) sprintf("Neoplastic vs non-neoplastic: chi-square %s; GEE (cluster-adjusted) %s",
nhan_p(p31_chi), nhan_p(p31_gee))
else sprintf("Tân sinh so với không tân sinh: χ² %s; GEE (hiệu chỉnh cụm) %s",
nhan_p(p31_chi), nhan_p(p31_gee))) +
theme_bao()
}
g31_vi <- ve_31("vi"); g31_en <- ve_31("en")
luu_hinh(g31_vi, "Bieu_do_3_1_Vi_tri_VN", 17.5, 10)
luu_hinh(g31_en, "Figure_3_1_Location_EN", 17.5, 10)
g31_vi
PHẦN 4. KẾT QUẢ MÔ BỆNH HỌC (STARD mục 21)
Mô tả phân bố chẩn đoán mô bệnh học (tiêu chuẩn tham chiếu) và tỷ lệ polyp tân sinh kèm KTC 95% Wilson.
tb33 <- data.frame(`Chẩn đoán mô bệnh học` = lab_mbh, `n (%)` = np(table(d$mbh_f)), check.names = FALSE)
tb33 <- rbind(tb33, data.frame(
`Chẩn đoán mô bệnh học` = "Tỷ lệ polyp tân sinh (u tuyến + SSL + carcinoma), % (KTC 95% Wilson)",
`n (%)` = fmt_ci(wilson(sum(d$ref_ts == 1), nrow(d))), check.names = FALSE))
dinh_dang_bang(flextable(tb33) |> bold(i = nrow(tb33)) |> align(j = 2, align = "center", part = "all"),
"Bảng 3.3. Đặc điểm mô bệnh học của polyp ≤ 5 mm",
sprintf("n = %d polyp.", nrow(d)))Chẩn đoán mô bệnh học | n (%) |
|---|---|
U tuyến | 129 (46.9) |
SSL | 1 (0.4) |
Tăng sản | 75 (27.3) |
Carcinoma | 0 (0.0) |
Viêm/khác | 70 (25.5) |
Tỷ lệ polyp tân sinh (u tuyến + SSL + carcinoma), % (KTC 95% Wilson) | 47.3 (41.5–53.2) |
n = 275 polyp. | |
PHẦN 5. ĐỘ CHÍNH XÁC CHẨN ĐOÁN CỦA BÁC SĨ NỘI SOI VÀ CADx (STARD mục 23–24)
5.1. Bảng 2 × 2 so với mô bệnh học
bang_2x2 <- function(test, ref, ten_test, ten_ref = "MBH") {
k <- dem(test, ref)
df <- data.frame(c(paste(ten_test, "– Tân sinh"), paste(ten_test, "– Không tân sinh"), "Tổng"),
c(k[["TP"]], k[["FN"]], k[["TP"]] + k[["FN"]]),
c(k[["FP"]], k[["TN"]], k[["FP"]] + k[["TN"]]),
c(k[["TP"]] + k[["FP"]], k[["FN"]] + k[["TN"]], sum(k)))
names(df) <- c(" ", paste0(ten_ref, ": Tân sinh"), paste0(ten_ref, ": Không tân sinh"), "Tổng"); df
}
dinh_dang_bang(flextable(bang_2x2(d_pair$bs_ts, d_pair$ref_ts, "Bác sĩ")) |> align(j = 2:4, align = "center", part = "all"),
"Bảng 3.4. Phân loại polyp giữa bác sĩ nội soi và mô bệnh học")
| MBH: Tân sinh | MBH: Không tân sinh | Tổng |
|---|---|---|---|
Bác sĩ – Tân sinh | 123 | 81 | 204 |
Bác sĩ – Không tân sinh | 7 | 64 | 71 |
Tổng | 130 | 145 | 275 |
dinh_dang_bang(flextable(bang_2x2(d_pair$cadx_ts, d_pair$ref_ts, "CADx")) |> align(j = 2:4, align = "center", part = "all"),
"Bảng 3.5. Phân loại polyp giữa CADx và mô bệnh học")
| MBH: Tân sinh | MBH: Không tân sinh | Tổng |
|---|---|---|---|
CADx – Tân sinh | 114 | 55 | 169 |
CADx – Không tân sinh | 16 | 90 | 106 |
Tổng | 130 | 145 | 275 |
5.2. Bảng 3.6 – So sánh giá trị chẩn đoán giữa bác sĩ và CADx (có p-value)
KTC 95%: bootstrap cụm theo bệnh nhân (B = 2000, seed = 2026).
p-value – hai cột:
Chỉ số p hiệu chỉnh cụm (kết quả chính) p không hiệu chỉnh (tham khảo, gói DTComPair)
Độ nhạy, độ đặc hiệu, độ chính xác GEE (cụm = bệnh nhân) McNemar chính xác
PPV, NPV Bootstrap cụm Điểm tổng quát của Leisenring, Alonzo & Pepe (2000)
LR+, LR− Bootstrap cụm Mô hình hồi quy của Gu & Pepe (2009)
tab.paired() của DTComPair tạo bảng bắt cặp (mô bệnh học × bác sĩ × CADx) từ ba cột 0/1.
# KTC bootstrap cụm (thứ tự tính giữ như phiên trước → KTC trùng khớp)
stat_cap <- function(x) {
a <- chi_so(x$bs_ts, x$ref_ts); b <- chi_so(x$cadx_ts, x$ref_ts)
c(setNames(a, paste0("BS_", names(a))), setNames(b, paste0("CADx_", names(b))),
setNames(b - a, paste0("D_", names(a))))
}
diem <- stat_cap(d_pair); bt <- boot_cum(d_pair, stat_cap, B, SEED); ci <- ci_q(bt)
# p hiệu chỉnh cụm
ts <- d_pair$ref_ts == 1; ks <- d_pair$ref_ts == 0
p_chinh <- c(Se = gee_p(d_pair[ts, ]), Sp = gee_p(d_pair[ks, ]), Acc = gee_p(d_pair),
PPV = p_boot(bt[, "D_PPV"]), NPV = p_boot(bt[, "D_NPV"]),
LRp = p_boot(bt[, "D_LRp"]), LRn = p_boot(bt[, "D_LRn"]))
# p không hiệu chỉnh (DTComPair + McNemar chính xác cho độ chính xác)
tab_cap <- tab.paired(d = d_pair$ref_ts, y1 = d_pair$bs_ts, y2 = d_pair$cadx_ts)
kq_sesp <- sesp.exactbinom(tab_cap)
kq_pv <- pv.gs(tab_cap)
kq_lr <- dlr.regtest(tab_cap)
p_phu <- c(Se = kq_sesp$sensitivity[["p.value"]], Sp = kq_sesp$specificity[["p.value"]],
Acc = mcnemar_cx(d_pair$bs_ts == d_pair$ref_ts, d_pair$cadx_ts == d_pair$ref_ts),
PPV = kq_pv$ppv[["p.value"]], NPV = kq_pv$npv[["p.value"]],
LRp = kq_lr$pdlr$p.value, LRn = kq_lr$ndlr$p.value)
ten_cs <- c(Se = "Độ nhạy", Sp = "Độ đặc hiệu", PPV = "Giá trị tiên đoán dương (PPV)",
NPV = "Giá trị tiên đoán âm (NPV)", Acc = "Độ chính xác",
LRp = "Tỷ số khả dĩ dương (LR+)", LRn = "Tỷ số khả dĩ âm (LR−)")
dong36 <- function(m) {
lr <- m %in% c("LRp", "LRn")
f <- function(pre) { v <- c(diem[paste0(pre, m)], ci[, paste0(pre, m)])
if (lr) sprintf("%.2f (%.2f–%.2f)", v[1], v[2], v[3]) else sprintf("%.1f (%.1f–%.1f)", 100 * v[1], 100 * v[2], 100 * v[3]) }
data.frame(`Chỉ số` = ten_cs[[m]], `Bác sĩ nội soi` = f("BS_"), `CADx` = f("CADx_"),
`Hiệu số CADx − BS` = f("D_"), `p hiệu chỉnh cụm` = fmtp(p_chinh[[m]]),
`p không hiệu chỉnh` = fmtp(p_phu[[m]]), check.names = FALSE)
}
tb36 <- bind_rows(lapply(names(ten_cs), dong36))
dinh_dang_bang(flextable(tb36) |> align(j = 2:6, align = "center", part = "all"),
"Bảng 3.6. So sánh khả năng phân loại polyp tân sinh và không tân sinh giữa bác sĩ nội soi và CADx",
c(sprintf("%% hoặc tỷ số (KTC 95%% bootstrap cụm theo bệnh nhân, B = %d); n = %d polyp/%d bệnh nhân.",
B, nrow(d_pair), n_distinct(d_pair$id_bn)),
"p hiệu chỉnh cụm: GEE logit, tương quan trao đổi (độ nhạy, độ đặc hiệu, độ chính xác); bootstrap cụm (PPV, NPV, LR+, LR−).",
"p không hiệu chỉnh: McNemar chính xác (độ nhạy, độ đặc hiệu, độ chính xác); điểm tổng quát Leisenring–Alonzo–Pepe (PPV, NPV); hồi quy Gu–Pepe (LR+, LR−)."))Chỉ số | Bác sĩ nội soi | CADx | Hiệu số CADx − BS | p hiệu chỉnh cụm | p không hiệu chỉnh |
|---|---|---|---|---|---|
Độ nhạy | 94.6 (90.1–98.0) | 87.7 (82.0–93.0) | -6.9 (-11.8–-2.2) | 0.011 | 0.022 |
Độ đặc hiệu | 44.1 (34.3–53.9) | 62.1 (53.9–70.0) | 17.9 (10.6–25.3) | <0.001 | <0.001 |
Giá trị tiên đoán dương (PPV) | 60.3 (52.7–68.0) | 67.5 (59.9–74.7) | 7.2 (3.4–11.0) | <0.001 | <0.001 |
Giá trị tiên đoán âm (NPV) | 90.1 (82.4–96.3) | 84.9 (77.8–91.4) | -5.2 (-11.5–1.4) | 0.106 | 0.131 |
Độ chính xác | 68.0 (61.6–74.0) | 74.2 (68.7–79.2) | 6.2 (1.5–10.9) | 0.008 | 0.019 |
Tỷ số khả dĩ dương (LR+) | 1.69 (1.43–2.06) | 2.31 (1.88–2.93) | 0.62 (0.28–1.09) | <0.001 | <0.001 |
Tỷ số khả dĩ âm (LR−) | 0.12 (0.04–0.23) | 0.20 (0.11–0.30) | 0.08 (-0.02–0.16) | 0.105 | 0.164 |
% hoặc tỷ số (KTC 95% bootstrap cụm theo bệnh nhân, B = 2000); n = 275 polyp/160 bệnh nhân. | |||||
p hiệu chỉnh cụm: GEE logit, tương quan trao đổi (độ nhạy, độ đặc hiệu, độ chính xác); bootstrap cụm (PPV, NPV, LR+, LR−). | |||||
p không hiệu chỉnh: McNemar chính xác (độ nhạy, độ đặc hiệu, độ chính xác); điểm tổng quát Leisenring–Alonzo–Pepe (PPV, NPV); hồi quy Gu–Pepe (LR+, LR−). | |||||
wil <- function(test, ref) { k <- dem(test, ref)
rbind(Se = wilson(k[["TP"]], k[["TP"]] + k[["FN"]]), Sp = wilson(k[["TN"]], k[["TN"]] + k[["FP"]]),
PPV = wilson(k[["TP"]], k[["TP"]] + k[["FP"]]), NPV = wilson(k[["TN"]], k[["TN"]] + k[["FN"]]),
Acc = wilson(k[["TP"]] + k[["TN"]], sum(k))) }
wbs <- wil(d_pair$bs_ts, d_pair$ref_ts); wcx <- wil(d_pair$cadx_ts, d_pair$ref_ts)
tbw <- data.frame(`Chỉ số` = ten_cs[rownames(wbs)], `Bác sĩ – Wilson` = apply(wbs, 1, fmt_ci),
`CADx – Wilson` = apply(wcx, 1, fmt_ci), check.names = FALSE)
dinh_dang_bang(flextable(tbw) |> align(j = 2:3, align = "center", part = "all"),
"Bảng 3.6b. Khoảng tin cậy Wilson (không hiệu chỉnh cụm) – tham khảo")Chỉ số | Bác sĩ – Wilson | CADx – Wilson |
|---|---|---|
Độ nhạy | 94.6 (89.3–97.4) | 87.7 (80.9–92.3) |
Độ đặc hiệu | 44.1 (36.3–52.3) | 62.1 (54.0–69.6) |
Giá trị tiên đoán dương (PPV) | 60.3 (53.4–66.8) | 67.5 (60.1–74.1) |
Giá trị tiên đoán âm (NPV) | 90.1 (81.0–95.1) | 84.9 (76.9–90.5) |
Độ chính xác | 68.0 (62.3–73.2) | 74.2 (68.7–79.0) |
# Xem đầy đủ kết quả DTComPair (để học và đối chiếu)
tab_cap## Two binary diagnostic tests (paired design)
##
## Test1: 'd_pair$bs_ts'
## Test2: 'd_pair$cadx_ts'
##
## Diseased:
## Test1 pos. Test1 neg. Total
## Test2 pos. 112 2 114
## Test2 neg. 11 5 16
## Total 123 7 130
##
## Non-diseased:
## Test1 pos. Test1 neg. Total
## Test2 pos. 51 4 55
## Test2 neg. 30 60 90
## Total 81 64 145kq_sesp## $sensitivity
## test1 test2 diff p.value
## 0.94615385 0.87692308 -0.06923077 0.02246094
##
## $specificity
## test1 test2 diff p.value
## 4.413793e-01 6.206897e-01 1.793103e-01 6.164890e-06
##
## $method
## [1] "exactbinom"kq_lr## $pdlr
## $pdlr$test1
## [1] 1.693732
##
## $pdlr$test2
## [1] 2.311888
##
## $pdlr$ratio
## [1] 0.7326186
##
## $pdlr$se.log
## [1] 0.0925148
##
## $pdlr$test.statistic
## [1] -3.36303
##
## $pdlr$p.value
## [1] 0.0007709198
##
## $pdlr$lcl
## [1] 0.6111238
##
## $pdlr$ucl
## [1] 0.8782672
##
##
## $ndlr
## $ndlr$test1
## [1] 0.1219952
##
## $ndlr$test2
## [1] 0.1982906
##
## $ndlr$ratio
## [1] 0.6152344
##
## $ndlr$se.log
## [1] 0.3492481
##
## $ndlr$test.statistic
## [1] -1.390851
##
## $ndlr$p.value
## [1] 0.1642706
##
## $ndlr$lcl
## [1] 0.3102845
##
## $ndlr$ucl
## [1] 1.219891
##
##
## $alpha
## [1] 0.05
##
## $method
## [1] "DLR regression model (regtest)"5.3. Biểu đồ 3.2 – Phân bố chẩn đoán của bác sĩ, CADx và mô bệnh học
Ngoặc trên từng nhóm: McNemar chính xác so sánh tỷ lệ bác sĩ và CADx chẩn đoán vào nhóm đó (vẽ bằng ggpubr::stat_pvalue_manual). Phụ đề: kiểm định Bowker cho toàn bộ bảng 3 × 3 (bác sĩ vs CADx; mỗi phương pháp vs mô bệnh học).
phan_bo <- bind_rows(data.frame(nguon = "MBH", lop = d_pair$ref3),
data.frame(nguon = "BS", lop = d_pair$bs3),
data.frame(nguon = "CADx", lop = d_pair$cadx3)) |>
count(nguon, lop, .drop = FALSE) |> group_by(nguon) |> mutate(tl = n / sum(n)) |> ungroup() |>
mutate(nguon = factor(nguon, levels = c("MBH", "BS", "CADx")))
# p McNemar cho từng nhóm (BS vs CADx)
p32_lop <- sapply(lv3, function(cl) mcnemar_cx(d_pair$bs3 == cl, d_pair$cadx3 == cl))
ngoac32 <- data.frame(lop = lv3, i = seq_along(lv3)) |>
mutate(group1 = "BS", group2 = "CADx", xmin = i, xmax = i + 0.8 / 3,
y = sapply(lv3, function(cl) max(phan_bo$tl[phan_bo$lop == cl & phan_bo$nguon != "MBH"])) + 0.09,
nhan = nhan_p(p32_lop))
bw_bc <- bowker(d_pair$bs3, d_pair$cadx3, lv3)[["p"]]
bw_bm <- bowker(d_pair$bs3, d_pair$ref3, lv3)[["p"]]
bw_cm <- bowker(d_pair$cadx3, d_pair$ref3, lv3)[["p"]]
ve_32 <- function(lang = "vi") {
en <- lang == "en"
nhan_lop <- if (en) c("Adenoma", "SSL", "Hyperplastic/non-neoplastic") else lv3
nhan_ng <- if (en) c(MBH = "Histopathology", BS = "Endoscopist", CADx = "CADx")
else c(MBH = "Mô bệnh học", BS = "Bác sĩ nội soi", CADx = "CADx")
ggplot(phan_bo, aes(lop, tl, fill = nguon)) +
geom_col(position = position_dodge(0.8), width = 0.75, colour = "white", linewidth = 0.2) +
geom_text(aes(label = sprintf("%d (%.1f%%)", n, 100 * tl)), position = position_dodge(0.8),
vjust = -0.4, size = 2.3, family = FONT) +
ggpubr::stat_pvalue_manual(ngoac32, label = "nhan", xmin = "xmin", xmax = "xmax",
y.position = "y", tip.length = 0.01, size = 2.4,
bracket.size = 0.3, inherit.aes = FALSE) +
scale_fill_manual(values = MAU_PP, labels = nhan_ng) +
scale_x_discrete(labels = nhan_lop) +
scale_y_continuous(labels = percent, limits = c(0, 1.05), expand = expansion(mult = c(0, 0.02))) +
labs(x = NULL, y = if (en) "Proportion of polyps" else "Tỷ lệ polyp",
subtitle = if (en) sprintf("Bowker test: endoscopist vs CADx %s;\nendoscopist vs histopathology %s; CADx vs histopathology %s",
nhan_p(bw_bc), nhan_p(bw_bm), nhan_p(bw_cm))
else sprintf("Kiểm định Bowker: bác sĩ vs CADx %s;\nbác sĩ vs mô bệnh học %s; CADx vs mô bệnh học %s",
nhan_p(bw_bc), nhan_p(bw_bm), nhan_p(bw_cm)),
caption = if (en) "Brackets: exact McNemar test, endoscopist vs CADx." else "Ngoặc: kiểm định McNemar chính xác, bác sĩ so với CADx.") +
theme_bao()
}
g32_vi <- ve_32("vi"); g32_en <- ve_32("en")
luu_hinh(g32_vi, "Bieu_do_3_2_Phan_bo_chan_doan_VN", 17.5, 10)
luu_hinh(g32_en, "Figure_3_2_Diagnosis_distribution_EN", 17.5, 10)
g32_vi
g32_en
5.4. Biểu đồ 3.3 – Chỉ số chẩn đoán của bác sĩ và CADx
Cột = ước lượng, thanh sai số = KTC 95% bootstrap cụm. Ngoặc p = p hiệu chỉnh cụm (kết quả chính của Bảng 3.6), vẽ bằng ggpubr::stat_pvalue_manual.
cs5 <- c("Se", "Sp", "PPV", "NPV", "Acc")
fp <- bind_rows(lapply(cs5, function(m) data.frame(
cs = m, pp = c("BS", "CADx"), est = diem[c(paste0("BS_", m), paste0("CADx_", m))],
lo = ci[1, c(paste0("BS_", m), paste0("CADx_", m))], hi = ci[2, c(paste0("BS_", m), paste0("CADx_", m))]))) |>
mutate(cs = factor(cs, levels = cs5), pp = factor(pp, levels = c("BS", "CADx")))
ngoac33 <- data.frame(cs = cs5, i = seq_along(cs5)) |>
mutate(group1 = "BS", group2 = "CADx", xmin = i - 0.2, xmax = i + 0.2,
y.position = sapply(cs5, function(m) max(fp$hi[fp$cs == m])) + 0.05,
p.lab = nhan_p(p_chinh[cs5]))
ve_33 <- function(lang = "vi") {
en <- lang == "en"
nhan_cs <- if (en) c("Sensitivity", "Specificity", "PPV", "NPV", "Accuracy")
else c("Độ nhạy", "Độ đặc hiệu", "PPV", "NPV", "Độ chính xác")
nhan_pp <- if (en) c(BS = "Endoscopist", CADx = "CADx") else c(BS = "Bác sĩ nội soi", CADx = "CADx")
ggplot(fp, aes(cs, est, fill = pp)) +
geom_col(position = position_dodge(0.8), width = 0.75, colour = "white", linewidth = 0.2) +
geom_errorbar(aes(ymin = lo, ymax = hi), position = position_dodge(0.8), width = 0.18, linewidth = 0.35) +
geom_text(aes(y = 0.03, label = sprintf("%.1f", 100 * est)), position = position_dodge(0.8),
vjust = 0, size = 2.4, colour = "white", fontface = "bold", family = FONT) +
ggpubr::stat_pvalue_manual(ngoac33, label = "p.lab", xmin = "xmin", xmax = "xmax",
y.position = "y.position", tip.length = 0.01, size = 2.4,
bracket.size = 0.3, inherit.aes = FALSE) +
geom_hline(yintercept = 0.9, linetype = "dashed", colour = "grey45", linewidth = 0.3) +
scale_fill_manual(values = MAU_PP[c("BS", "CADx")], labels = nhan_pp) +
scale_x_discrete(labels = nhan_cs) +
scale_y_continuous(labels = percent, limits = c(0, 1.12), breaks = seq(0, 1, 0.2),
expand = expansion(mult = c(0, 0))) +
labs(x = NULL, y = if (en) "Estimate (95% CI)" else "Ước lượng (KTC 95%)",
caption = if (en) sprintf("Error bars: 95%% CI, patient-level cluster bootstrap (B = %d). p: GEE (sensitivity, specificity, accuracy) or\ncluster bootstrap (PPV, NPV). Dashed line: 90%% (PIVI reference).", B)
else sprintf("Thanh sai số: KTC 95%% bootstrap cụm theo bệnh nhân (B = %d). p: GEE (độ nhạy, độ đặc hiệu, độ chính xác)\nhoặc bootstrap cụm (PPV, NPV). Đường đứt đoạn: mốc 90%% (tham khảo PIVI).", B)) +
theme_bao()
}
g33_vi <- ve_33("vi"); g33_en <- ve_33("en")
luu_hinh(g33_vi, "Bieu_do_3_3_Chi_so_chan_doan_VN", 17.5, 10)
luu_hinh(g33_en, "Figure_3_3_Diagnostic_performance_EN", 17.5, 10)
g33_vi
g33_en
PHẦN 6. MỨC ĐỘ ĐỒNG THUẬN GIỮA BÁC SĨ NỘI SOI VÀ CADx
Ba loại p-value được báo cáo:
Kappa khác 0? – kiểm định z (phương sai dưới H0, Fleiss–Cohen–Everitt 1969).
Bác sĩ và CADx có xu hướng chẩn đoán “tân sinh” khác nhau? – McNemar chính xác (2 nhóm) và Bowker (3 nhóm).
Kappa của CADx với mô bệnh học có khác Kappa của bác sĩ với mô bệnh học? – hiệu số hai Kappa, KTC 95% và p bằng bootstrap cụm theo bệnh nhân.
Phân mức Kappa theo Landis & Koch (1977).
dinh_dang_bang(flextable(bang_2x2(d_pair$cadx_ts, d_pair$bs_ts, "CADx", "BS")) |>
align(j = 2:4, align = "center", part = "all"),
"Bảng 3.7. Kết quả phân loại polyp tân sinh và không tân sinh giữa bác sĩ nội soi và CADx")
| BS: Tân sinh | BS: Không tân sinh | Tổng |
|---|---|---|---|
CADx – Tân sinh | 163 | 6 | 169 |
CADx – Không tân sinh | 41 | 65 | 106 |
Tổng | 204 | 71 | 275 |
stat_k <- function(x) {
k2 <- kappa_z(x$bs_ts, x$cadx_ts, 0:1); k3 <- kappa_z(x$bs3, x$cadx3, lv3)
kb <- kappa_z(x$bs_ts, x$ref_ts, 0:1)[["k"]]; kc <- kappa_z(x$cadx_ts, x$ref_ts, 0:1)[["k"]]
c(po = k2[["po"]], k = k2[["k"]], po3 = k3[["po"]], k3 = k3[["k"]],
k_bs_ref = kb, k_cx_ref = kc, d_k = kc - kb)
}
kd <- stat_k(d_pair); bk <- boot_cum(d_pair, stat_k, B, SEED + 1); kb <- ci_q(bk)
# p-value
kz2 <- kappa_z(d_pair$bs_ts, d_pair$cadx_ts, 0:1)
kz3 <- kappa_z(d_pair$bs3, d_pair$cadx3, lv3)
kzb <- kappa_z(d_pair$bs_ts, d_pair$ref_ts, 0:1)
kzc <- kappa_z(d_pair$cadx_ts, d_pair$ref_ts, 0:1)
p_mc2 <- mcnemar_cx(d_pair$bs_ts == 1, d_pair$cadx_ts == 1)
p_bw3 <- bowker(d_pair$bs3, d_pair$cadx3, lv3)[["p"]]
p_dk <- p_boot(bk[, "d_k"])
muc_k <- function(k) as.character(cut(k, c(-Inf, 0, .2, .4, .6, .8, 1.01),
labels = c("Kém", "Rất yếu", "Yếu", "Trung bình", "Khá", "Rất tốt")))
fk <- function(nm, pct = FALSE) if (pct) sprintf("%.1f%% (%.1f–%.1f)", 100 * kd[nm], 100 * kb[1, nm], 100 * kb[2, nm]) else
sprintf("%.2f (%.2f–%.2f)", kd[nm], kb[1, nm], kb[2, nm])
tl_bs <- mean(d_pair$bs_ts == 1); tl_cx <- mean(d_pair$cadx_ts == 1)
tb38 <- data.frame(
`Chỉ số đánh giá` = c("Tỷ lệ đồng thuận thực tế p₀ – 2 nhóm",
"Kappa bác sĩ – CADx (2 nhóm)",
sprintf("Tỷ lệ chẩn đoán tân sinh: bác sĩ %.1f%% so với CADx %.1f%%", 100 * tl_bs, 100 * tl_cx),
"Tỷ lệ đồng thuận thực tế p₀ – 3 nhóm",
"Kappa bác sĩ – CADx (3 nhóm)",
"Đối xứng của bảng bác sĩ × CADx (3 nhóm)",
"Kappa bác sĩ – mô bệnh học (2 nhóm)",
"Kappa CADx – mô bệnh học (2 nhóm)",
"Hiệu số Kappa (CADx – MBH) − (bác sĩ – MBH)"),
`Giá trị (KTC 95%)` = c(fk("po", TRUE), fk("k"), sprintf("%.1f%% / %.1f%%", 100 * tl_bs, 100 * tl_cx),
fk("po3", TRUE), fk("k3"), "—", fk("k_bs_ref"), fk("k_cx_ref"), fk("d_k")),
`Mức độ (Landis & Koch)` = c("", muc_k(kd["k"]), "", "", muc_k(kd["k3"]), "", muc_k(kd["k_bs_ref"]),
muc_k(kd["k_cx_ref"]), ""),
p = fmtp(c(NA, kz2[["p"]], p_mc2, NA, kz3[["p"]], p_bw3, kzb[["p"]], kzc[["p"]], p_dk)),
`Kiểm định` = c("", "z (Kappa = 0)", "McNemar chính xác", "", "z (Kappa = 0)", "Bowker",
"z (Kappa = 0)", "z (Kappa = 0)", "Bootstrap cụm"),
check.names = FALSE)
dinh_dang_bang(flextable(tb38) |> align(j = 2:5, align = "center", part = "all"),
"Bảng 3.8. Mức độ đồng thuận trong phân loại polyp giữa bác sĩ nội soi và CADx",
c(sprintf("KTC 95%% bootstrap cụm theo bệnh nhân (B = %d). Phân mức Kappa theo Landis & Koch (1977).", B),
"z (Kappa = 0): phương sai dưới giả thuyết không theo Fleiss, Cohen & Everitt (1969)."))Chỉ số đánh giá | Giá trị (KTC 95%) | Mức độ (Landis & Koch) | p | Kiểm định |
|---|---|---|---|---|
Tỷ lệ đồng thuận thực tế p₀ – 2 nhóm | 82.9% (78.6–87.3) | — | ||
Kappa bác sĩ – CADx (2 nhóm) | 0.62 (0.51–0.72) | Khá | <0.001 | z (Kappa = 0) |
Tỷ lệ chẩn đoán tân sinh: bác sĩ 74.2% so với CADx 61.5% | 74.2% / 61.5% | <0.001 | McNemar chính xác | |
Tỷ lệ đồng thuận thực tế p₀ – 3 nhóm | 82.2% (77.8–86.7) | — | ||
Kappa bác sĩ – CADx (3 nhóm) | 0.61 (0.51–0.71) | Khá | <0.001 | z (Kappa = 0) |
Đối xứng của bảng bác sĩ × CADx (3 nhóm) | — | <0.001 | Bowker | |
Kappa bác sĩ – mô bệnh học (2 nhóm) | 0.38 (0.27–0.49) | Yếu | <0.001 | z (Kappa = 0) |
Kappa CADx – mô bệnh học (2 nhóm) | 0.49 (0.39–0.59) | Trung bình | <0.001 | z (Kappa = 0) |
Hiệu số Kappa (CADx – MBH) − (bác sĩ – MBH) | 0.11 (0.03–0.19) | 0.013 | Bootstrap cụm | |
KTC 95% bootstrap cụm theo bệnh nhân (B = 2000). Phân mức Kappa theo Landis & Koch (1977). | ||||
z (Kappa = 0): phương sai dưới giả thuyết không theo Fleiss, Cohen & Everitt (1969). | ||||
**PHẦN 7. PHÂN TÍCH 3 NHÓM MÔ BỆNH HỌC*
Ba nhóm: U tuyến (gộp carcinoma), SSL, Tăng sản/không tân sinh (gộp viêm/khác). Độ nhạy, độ đặc hiệu từng nhóm tính theo kiểu “một nhóm so với hai nhóm còn lại”. p so sánh bác sĩ với CADx: McNemar chính xác. Nhóm có < 5 polyp (SSL) không ước lượng KTC và không kiểm định.
bang3 <- function(test, ten_test) {
t <- table(factor(test, lv3), d_pair$ref3); m <- as.data.frame.matrix(t)
m <- cbind(` ` = paste(ten_test, "–", rownames(m)), m, `Tổng` = rowSums(t))
rbind(m, c(" " = "Tổng", colSums(t), sum(t))) }
dinh_dang_bang(flextable(bang3(d_pair$bs3, "Bác sĩ")) |>
add_header_row(values = c("", "Mô bệnh học", ""), colwidths = c(1, 3, 1)) |>
align(j = 2:5, align = "center", part = "all"),
"Bảng 3.9. Phân loại 3 nhóm giữa bác sĩ nội soi và mô bệnh học")Mô bệnh học | ||||
|---|---|---|---|---|
| U tuyến | SSL | Tăng sản/không tân sinh | Tổng |
Bác sĩ – U tuyến | 122 | 0 | 79 | 201 |
Bác sĩ – SSL | 0 | 1 | 2 | 3 |
Bác sĩ – Tăng sản/không tân sinh | 7 | 0 | 64 | 71 |
Tổng | 129 | 1 | 145 | 275 |
dinh_dang_bang(flextable(bang3(d_pair$cadx3, "CADx")) |>
add_header_row(values = c("", "Mô bệnh học", ""), colwidths = c(1, 3, 1)) |>
align(j = 2:5, align = "center", part = "all"),
"Bảng 3.10. Phân loại 3 nhóm giữa CADx và mô bệnh học")Mô bệnh học | ||||
|---|---|---|---|---|
| U tuyến | SSL | Tăng sản/không tân sinh | Tổng |
CADx – U tuyến | 113 | 0 | 55 | 168 |
CADx – SSL | 0 | 1 | 0 | 1 |
CADx – Tăng sản/không tân sinh | 16 | 0 | 90 | 106 |
Tổng | 129 | 1 | 145 | 275 |
dinh_dang_bang(flextable(bang3(d_pair$cadx3, "CADx")) |>
add_header_row(values = c("", "Mô bệnh học", ""), colwidths = c(1, 3, 1)) |>
align(j = 2:5, align = "center", part = "all"),
"Bảng 3.10. Phân loại 3 nhóm giữa CADx và mô bệnh học")Mô bệnh học | ||||
|---|---|---|---|---|
| U tuyến | SSL | Tăng sản/không tân sinh | Tổng |
CADx – U tuyến | 113 | 0 | 55 | 168 |
CADx – SSL | 0 | 1 | 0 | 1 |
CADx – Tăng sản/không tân sinh | 16 | 0 | 90 | 106 |
Tổng | 129 | 1 | 145 | 275 |
stat3 <- function(x) { out <- c()
for (pp in c("bs3", "cadx3")) {
for (cl in lv3) { s <- chi_so(as.integer(x[[pp]] == cl), as.integer(x$ref3 == cl))
out[paste(pp, cl, "Se")] <- s[["Se"]]; out[paste(pp, cl, "Sp")] <- s[["Sp"]] }
out[paste(pp, "Acc3")] <- mean(x[[pp]] == x$ref3) }
out }
s3 <- stat3(d_pair); b3 <- ci_q(boot_cum(d_pair, stat3, B, SEED + 2))
n_lop <- table(d_pair$ref3)
f3 <- function(nm, cl = NULL) {
if (is.na(s3[nm])) return("—")
if (!is.null(cl) && n_lop[[cl]] < 5) return(sprintf("%.1f (không ước lượng KTC)", 100 * s3[nm]))
sprintf("%.1f (%.1f–%.1f)", 100 * s3[nm], 100 * b3[1, nm], 100 * b3[2, nm]) }
# p McNemar chính xác
p3 <- function(cl, loai) {
if (n_lop[[cl]] < 5) return(NA)
if (loai == "Se") { i <- d_pair$ref3 == cl; mcnemar_cx(d_pair$bs3[i] == cl, d_pair$cadx3[i] == cl) }
else { i <- d_pair$ref3 != cl; mcnemar_cx(d_pair$bs3[i] != cl, d_pair$cadx3[i] != cl) }
}
tb311 <- bind_rows(lapply(lv3, function(cl) data.frame(
`Nhóm mô bệnh học` = sprintf("%s (n = %d)", cl, n_lop[[cl]]),
`Độ nhạy – BS` = f3(paste("bs3", cl, "Se"), cl), `Độ nhạy – CADx` = f3(paste("cadx3", cl, "Se"), cl),
`p (độ nhạy)` = fmtp(p3(cl, "Se")),
`Độ đặc hiệu – BS` = f3(paste("bs3", cl, "Sp"), cl), `Độ đặc hiệu – CADx` = f3(paste("cadx3", cl, "Sp"), cl),
`p (độ đặc hiệu)` = fmtp(p3(cl, "Sp")), check.names = FALSE)))
p_acc3 <- mcnemar_cx(d_pair$bs3 == d_pair$ref3, d_pair$cadx3 == d_pair$ref3)
tb311 <- rbind(tb311, data.frame(
`Nhóm mô bệnh học` = "Độ chính xác chung 3 nhóm", `Độ nhạy – BS` = f3("bs3 Acc3"),
`Độ nhạy – CADx` = f3("cadx3 Acc3"), `p (độ nhạy)` = fmtp(p_acc3),
`Độ đặc hiệu – BS` = "", `Độ đặc hiệu – CADx` = "", `p (độ đặc hiệu)` = "", check.names = FALSE))
chu_thich3 <- c("% (KTC 95% bootstrap cụm). Mỗi nhóm so với hai nhóm còn lại. Carcinoma gộp vào u tuyến; viêm/khác gộp vào không tân sinh.",
"p: kiểm định McNemar chính xác, bác sĩ so với CADx (không hiệu chỉnh cụm). Dòng độ chính xác chung: giá trị nằm ở cột độ nhạy.")
if (n_lop[["SSL"]] < 10) chu_thich3 <- c(chu_thich3,
sprintf("Nhóm SSL chỉ có %d polyp – không ước lượng KTC, không kiểm định; số liệu chỉ mang tính mô tả.", n_lop[["SSL"]]))
dinh_dang_bang(flextable(tb311) |> bold(i = nrow(tb311)) |> align(j = 2:7, align = "center", part = "all"),
"Bảng 3.11. Độ nhạy, độ đặc hiệu theo từng nhóm mô bệnh học và độ chính xác 3 nhóm", chu_thich3)Nhóm mô bệnh học | Độ nhạy – BS | Độ nhạy – CADx | p (độ nhạy) | Độ đặc hiệu – BS | Độ đặc hiệu – CADx | p (độ đặc hiệu) |
|---|---|---|---|---|---|---|
U tuyến (n = 129) | 94.6 (90.0–98.1) | 87.6 (81.7–92.5) | 0.022 | 45.9 (36.6–55.7) | 62.3 (53.9–70.4) | <0.001 |
SSL (n = 1) | 100.0 (không ước lượng KTC) | 100.0 (không ước lượng KTC) | — | 99.3 (không ước lượng KTC) | 100.0 (không ước lượng KTC) | — |
Tăng sản/không tân sinh (n = 145) | 44.1 (34.8–53.7) | 62.1 (53.9–70.1) | <0.001 | 94.6 (90.1–98.1) | 87.7 (81.9–92.5) | 0.022 |
Độ chính xác chung 3 nhóm | 68.0 (61.4–74.4) | 74.2 (68.6–79.1) | 0.019 | |||
% (KTC 95% bootstrap cụm). Mỗi nhóm so với hai nhóm còn lại. Carcinoma gộp vào u tuyến; viêm/khác gộp vào không tân sinh. | ||||||
p: kiểm định McNemar chính xác, bác sĩ so với CADx (không hiệu chỉnh cụm). Dòng độ chính xác chung: giá trị nằm ở cột độ nhạy. | ||||||
Nhóm SSL chỉ có 1 polyp – không ước lượng KTC, không kiểm định; số liệu chỉ mang tính mô tả. | ||||||
7.1. Biểu đồ 3.4 – Ma trận nhầm lẫn 3 nhóm
k3_bs <- kappa_z(d_pair$bs3, d_pair$ref3, lv3); k3_cx <- kappa_z(d_pair$cadx3, d_pair$ref3, lv3)
ve_34 <- function(lang = "vi") {
en <- lang == "en"
nhan_lop <- if (en) c("Adenoma", "SSL", "Hyperplastic/\nnon-neoplastic") else c("U tuyến", "SSL", "Tăng sản/\nkhông tân sinh")
hm <- function(test, tieu_de, truc_y) {
as.data.frame(table(Test = factor(test, lv3), MBH = d_pair$ref3)) |> group_by(MBH) |>
mutate(tl = if (sum(Freq) > 0) Freq / sum(Freq) else 0) |> ungroup() |> # % theo cột mô bệnh học
ggplot(aes(MBH, Test, fill = tl)) + geom_tile(colour = "white", linewidth = 0.6) +
geom_text(aes(label = sprintf("%d\n(%.1f%%)", Freq, 100 * tl), colour = tl > 0.55),
size = 2.5, lineheight = 0.9, family = FONT, show.legend = FALSE) +
scale_colour_manual(values = c(`TRUE` = "white", `FALSE` = "black")) +
scale_fill_gradient(low = "#F2F2F2", high = "#0072B2", limits = c(0, 1), labels = percent,
name = if (en) "Column %" else "% theo cột") +
scale_x_discrete(labels = setNames(nhan_lop, lv3)) +
scale_y_discrete(limits = rev(lv3), labels = setNames(nhan_lop, lv3)) +
labs(x = if (en) "Histopathology" else "Mô bệnh học", y = truc_y, subtitle = tieu_de) +
theme_bao() + theme(legend.position = "right", legend.title = element_text(size = 7),
legend.key.height = unit(6, "mm"), axis.line = element_blank(),
axis.ticks = element_blank(), plot.subtitle = element_text(face = "bold", size = 8))
}
t_bs <- if (en) sprintf("A. Endoscopist\nAccuracy %.1f%%; κ = %.2f", 100 * s3[["bs3 Acc3"]], k3_bs[["k"]])
else sprintf("A. Bác sĩ nội soi\nĐộ chính xác %.1f%%; κ = %.2f", 100 * s3[["bs3 Acc3"]], k3_bs[["k"]])
t_cx <- if (en) sprintf("B. CADx\nAccuracy %.1f%%; κ = %.2f", 100 * s3[["cadx3 Acc3"]], k3_cx[["k"]])
else sprintf("B. CADx\nĐộ chính xác %.1f%%; κ = %.2f", 100 * s3[["cadx3 Acc3"]], k3_cx[["k"]])
(hm(d_pair$bs3, t_bs, if (en) "Endoscopist" else "Bác sĩ nội soi") +
hm(d_pair$cadx3, t_cx, "CADx")) +
plot_layout(guides = "collect") &
theme(legend.position = "right")
}
g34_vi <- ve_34("vi") + plot_annotation(caption = sprintf("So sánh độ chính xác 3 nhóm, bác sĩ vs CADx: McNemar chính xác %s. κ: Kappa Cohen so với mô bệnh học (3 nhóm).", nhan_p(p_acc3)),
theme = theme(plot.caption = element_text(size = 7, family = FONT, hjust = 0)))
g34_en <- ve_34("en") + plot_annotation(caption = sprintf("Three-class accuracy, endoscopist vs CADx: exact McNemar %s. κ: Cohen kappa vs histopathology (3 classes).", nhan_p(p_acc3)),
theme = theme(plot.caption = element_text(size = 7, family = FONT, hjust = 0)))
luu_hinh(g34_vi, "Bieu_do_3_4_Ma_tran_nham_lan_VN", 17.5, 9)
luu_hinh(g34_en, "Figure_3_4_Confusion_matrix_EN", 17.5, 9)
g34_vi
g34_en
PHẦN 8. PHÂN TÍCH DƯỚI NHÓM THEO VỊ TRÍ (thăm dò)
So sánh vùng trực tràng – sigma với các đoạn còn lại (liên quan chiến lược “chẩn đoán và để lại” theo ngưỡng PIVI: NPV ≥ 90% cho polyp ≤ 5 mm vùng trực tràng – sigma). p so sánh bác sĩ với CADx trong từng vùng: McNemar chính xác (độ nhạy, độ đặc hiệu), Leisenring (DTComPair::pv.gs, NPV). Phân tích thăm dò, không hiệu chỉnh cụm.
tb312 <- bind_rows(lapply(levels(d_pair$vung), function(v) {
x <- d_pair[d_pair$vung == v, ]
kb_ <- dem(x$bs_ts, x$ref_ts); kc_ <- dem(x$cadx_ts, x$ref_ts)
pv_v <- pv.gs(tab.paired(d = x$ref_ts, y1 = x$bs_ts, y2 = x$cadx_ts))
t_ <- x$ref_ts == 1; k_ <- x$ref_ts == 0
data.frame(
`Vùng` = sprintf("%s\nn = %d (tân sinh %d/không %d)", v, nrow(x), sum(t_), sum(k_)),
`Chỉ số` = c("Độ nhạy", "Độ đặc hiệu", "NPV"),
`Bác sĩ nội soi` = c(fmt_ci(wilson(kb_[["TP"]], kb_[["TP"]] + kb_[["FN"]])),
fmt_ci(wilson(kb_[["TN"]], kb_[["TN"]] + kb_[["FP"]])),
fmt_ci(wilson(kb_[["TN"]], kb_[["TN"]] + kb_[["FN"]]))),
`CADx` = c(fmt_ci(wilson(kc_[["TP"]], kc_[["TP"]] + kc_[["FN"]])),
fmt_ci(wilson(kc_[["TN"]], kc_[["TN"]] + kc_[["FP"]])),
fmt_ci(wilson(kc_[["TN"]], kc_[["TN"]] + kc_[["FN"]]))),
p = fmtp(c(mcnemar_cx(x$bs_ts[t_] == 1, x$cadx_ts[t_] == 1),
mcnemar_cx(x$bs_ts[k_] == 0, x$cadx_ts[k_] == 0),
pv_v$npv[["p.value"]])),
`Kiểm định` = c("McNemar chính xác", "McNemar chính xác", "Leisenring"),
check.names = FALSE)
}))
dinh_dang_bang(flextable(tb312) |> merge_v(j = 1) |> valign(j = 1, valign = "top") |>
align(j = 3:6, align = "center", part = "all"),
"Bảng 3.12. Giá trị chẩn đoán theo vị trí polyp",
c("% (KTC 95% Wilson). Phân tích thăm dò, không hiệu chỉnh cụm.",
"NPV vùng trực tràng – sigma so sánh tham khảo với ngưỡng PIVI (≥ 90%); không ghi nhận mức độ tin cậy của chẩn đoán nên không đánh giá đầy đủ tiêu chuẩn PIVI."))Vùng | Chỉ số | Bác sĩ nội soi | CADx | p | Kiểm định |
|---|---|---|---|---|---|
Trực tràng – sigma | Độ nhạy | 88.6 (74.0–95.5) | 91.4 (77.6–97.0) | 1.000 | McNemar chính xác |
Độ đặc hiệu | 54.5 (43.5–65.2) | 72.7 (61.9–81.4) | <0.001 | McNemar chính xác | |
NPV | 91.3 (79.7–96.6) | 94.9 (86.1–98.3) | 0.284 | Leisenring | |
Đại tràng phải, ngang, trái | Độ nhạy | 96.8 (91.1–98.9) | 86.3 (78.0–91.8) | 0.002 | McNemar chính xác |
Độ đặc hiệu | 32.4 (22.4–44.2) | 50.0 (38.4–61.6) | 0.012 | McNemar chính xác | |
NPV | 88.0 (70.0–95.8) | 72.3 (58.2–83.1) | 0.018 | Leisenring | |
% (KTC 95% Wilson). Phân tích thăm dò, không hiệu chỉnh cụm. | |||||
NPV vùng trực tràng – sigma so sánh tham khảo với ngưỡng PIVI (≥ 90%); không ghi nhận mức độ tin cậy của chẩn đoán nên không đánh giá đầy đủ tiêu chuẩn PIVI. | |||||
8.1. Đối chiếu với ngưỡng PIVI (ASGE) và SODA (ESGE)
Ngưỡng tham chiếu (đã đối chiếu với tài liệu gốc):
Chiến lược Tổ chức Ngưỡng Phạm vi
Để lại tại chỗ (diagnose-and-leave) ASGE PIVI (Rex 2011) NPV ≥ 90% Polyp ≤ 5 mm trực tràng – sigma
Cắt và bỏ (resect-and-discard) ASGE PIVI (Rex 2011) Đồng thuận khoảng theo dõi ≥ 90% Toàn bộ đại tràng
Để lại tại chỗ ESGE SODA (Houwen 2022) Độ nhạy ≥ 90% và độ đặc hiệu ≥ 80% Polyp 1–5 mm trực tràng – sigma
Cắt và bỏ ESGE SODA (Houwen 2022) Độ nhạy ≥ 80% và độ đặc hiệu ≥ 80% Polyp 1–5 mm toàn bộ đại tràng
Quy tắc đánh giá (so với ngưỡng):
✓ Đạt: cận dưới KTC 95% ≥ ngưỡng.
≈ Đạt ước lượng điểm: ước lượng điểm ≥ ngưỡng nhưng cận dưới KTC 95% < ngưỡng.
✗ Chưa đạt: ước lượng điểm < ngưỡng. Với tiêu chí gồm 2 chỉ số (SODA), chỉ cần 1 chỉ số chưa đạt thì tiêu chí chưa đạt.
KTC 95% tính bằng bootstrap cụm theo bệnh nhân (B = 2000); tập trực tràng – sigma được bootstrap riêng. p so sánh bác sĩ với CADx: GEE (độ nhạy, độ đặc hiệu, độ chính xác), bootstrap cụm (PPV, NPV). Tiêu chí PIVI “cắt và bỏ” không đánh giá được vì cần khoảng theo dõi sau cắt polyp và mô bệnh học của polyp > 5 mm – không có trong bộ dữ liệu.
NGUONG_PIVI <- 0.90; NGUONG_SODA <- 0.80
# Tập trực tràng – sigma và bootstrap cụm riêng
d_rs <- d_pair |> filter(vung == "Trực tràng – sigma")
diem_rs <- stat_cap(d_rs); bt_rs <- boot_cum(d_rs, stat_cap, B, SEED + 3); ci_rs <- ci_q(bt_rs)
p_rs <- c(Se = gee_p(d_rs[d_rs$ref_ts == 1, ]), Sp = gee_p(d_rs[d_rs$ref_ts == 0, ]), Acc = gee_p(d_rs),
PPV = p_boot(bt_rs[, "D_PPV"]), NPV = p_boot(bt_rs[, "D_NPV"]))
# Số thập phân kiểu Việt Nam (dấu phẩy) cho bản tiếng Việt
so_vn <- function(x, k = 1) sub(".", ",", sprintf("%.*f", k, x), fixed = TRUE)
fmtp_vn <- function(p) ifelse(is.na(p), "—", ifelse(p < 0.001, "< 0,001", sub(".", ",", sprintf("%.3f", p), fixed = TRUE)))
ci_vn <- function(e, l, h) sprintf("%s (%s–%s)", so_vn(100 * e), so_vn(100 * l), so_vn(100 * h))
# Hàm đánh giá một chỉ số so với ngưỡng: 2 = đạt, 1 = đạt ước lượng điểm, 0 = chưa đạt
muc_dat <- function(est, lo, thr) ifelse(lo >= thr, 2L, ifelse(est >= thr, 1L, 0L))
nhan_dat <- function(muc, thieu = NULL) switch(as.character(muc),
`2` = "✓ Đạt",
`1` = "≈ Đạt ước lượng điểm\n(cận dưới KTC < ngưỡng)",
`0` = paste0("✗ Chưa đạt", if (length(thieu)) paste0(" (", paste(thieu, collapse = ", "), ")") else ""))
# Lấy ước lượng + KTC của 1 chỉ số cho 1 phương pháp
lay <- function(diem_x, ci_x, pp, m) c(est = diem_x[[paste0(pp, "_", m)]], lo = ci_x[1, paste0(pp, "_", m)],
hi = ci_x[2, paste0(pp, "_", m)])
# Đánh giá tiêu chí (1 hoặc 2 chỉ số) cho 1 phương pháp
danh_gia <- function(diem_x, ci_x, pp, tieu_chi) { # tieu_chi = list(Se = 0.9, Sp = 0.8)
muc <- sapply(names(tieu_chi), function(m) { v <- lay(diem_x, ci_x, pp, m); muc_dat(v[["est"]], v[["lo"]], tieu_chi[[m]]) })
ten_vn <- c(Se = "độ nhạy", Sp = "độ đặc hiệu", NPV = "NPV")
list(muc = min(muc), nhan = nhan_dat(min(muc), if (min(muc) == 0) ten_vn[names(muc)[muc == 0]] else NULL))
}
gia_tri <- function(diem_x, ci_x, pp, ms) paste(sapply(ms, function(m) { v <- lay(diem_x, ci_x, pp, m)
paste0(c(Se = "Se ", Sp = "Sp ", NPV = "")[[m]], ci_vn(v[["est"]], v[["lo"]], v[["hi"]])) }), collapse = "\n")
n_rs <- sprintf("n = %d polyp/%d bệnh nhân", nrow(d_rs), n_distinct(d_rs$id_bn))
n_all <- sprintf("n = %d polyp/%d bệnh nhân", nrow(d_pair), n_distinct(d_pair$id_bn))
tieu_chi_ds <- list(
list(ten = "PIVI – để lại tại chỗ: NPV ≥ 90%", pv = paste("Trực tràng – sigma,", n_rs),
dx = diem_rs, cx = ci_rs, tc = list(NPV = NGUONG_PIVI)),
list(ten = "SODA – để lại tại chỗ: Se ≥ 90% và Sp ≥ 80%", pv = paste("Trực tràng – sigma,", n_rs),
dx = diem_rs, cx = ci_rs, tc = list(Se = NGUONG_PIVI, Sp = NGUONG_SODA)),
list(ten = "SODA – cắt và bỏ: Se ≥ 80% và Sp ≥ 80%", pv = paste("Toàn bộ đại tràng,", n_all),
dx = diem, cx = ci, tc = list(Se = NGUONG_SODA, Sp = NGUONG_SODA)))
tb312b <- bind_rows(lapply(tieu_chi_ds, function(t) {
a <- danh_gia(t$dx, t$cx, "BS", t$tc); b <- danh_gia(t$dx, t$cx, "CADx", t$tc)
data.frame(`Tiêu chí` = paste0(t$ten, "\n", t$pv),
`Bác sĩ nội soi` = gia_tri(t$dx, t$cx, "BS", names(t$tc)),
`CADx` = gia_tri(t$dx, t$cx, "CADx", names(t$tc)),
`Đánh giá – bác sĩ` = a$nhan, `Đánh giá – CADx` = b$nhan,
muc_bs = a$muc, muc_cx = b$muc, check.names = FALSE) }))
tb312b <- bind_rows(tb312b, data.frame(
`Tiêu chí` = "PIVI – cắt và bỏ: đồng thuận khoảng theo dõi ≥ 90%\nToàn bộ đại tràng",
`Bác sĩ nội soi` = "—", `CADx` = "—",
`Đánh giá – bác sĩ` = "Không đánh giá được", `Đánh giá – CADx` = "Không đánh giá được",
muc_bs = NA, muc_cx = NA, check.names = FALSE))mau_dat <- c(`2` = "#E2F0D9", `1` = "#FFF2CC", `0` = "#F8D7DA")
ft312b <- flextable(tb312b, col_keys = c("Tiêu chí", "Bác sĩ nội soi", "CADx", "Đánh giá – bác sĩ", "Đánh giá – CADx")) |>
align(j = 2:5, align = "center", part = "all")
for (i in seq_len(nrow(tb312b))) {
if (!is.na(tb312b$muc_bs[i])) ft312b <- bg(ft312b, i = i, j = "Đánh giá – bác sĩ", bg = mau_dat[[as.character(tb312b$muc_bs[i])]])
if (!is.na(tb312b$muc_cx[i])) ft312b <- bg(ft312b, i = i, j = "Đánh giá – CADx", bg = mau_dat[[as.character(tb312b$muc_cx[i])]])
}
ft312b <- dinh_dang_bang(ft312b |> bold(j = 4:5, part = "body"),
"Bảng 3.12b. Đối chiếu giá trị chẩn đoán của bác sĩ nội soi và CADx với ngưỡng PIVI (ASGE) và SODA (ESGE)",
c("% (KTC 95% bootstrap cụm theo bệnh nhân, B = 2000). Se: độ nhạy; Sp: độ đặc hiệu; NPV: giá trị tiên đoán âm – tính cho polyp tân sinh (u tuyến + SSL).",
"✓ Đạt: cận dưới KTC 95% ≥ ngưỡng; ≈ Đạt ước lượng điểm: ước lượng ≥ ngưỡng nhưng cận dưới KTC < ngưỡng; ✗ Chưa đạt: ước lượng < ngưỡng (ghi chỉ số chưa đạt).",
"Ngưỡng: ASGE PIVI (Rex và cs., Gastrointest Endosc 2011;73:419–22); ESGE SODA (Houwen và cs., Endoscopy 2022;54:88–99). Các ngưỡng quy định cho chẩn đoán với mức độ tin cậy cao; nghiên cứu không ghi nhận mức độ tin cậy nên đây là đối chiếu tham khảo.",
"PIVI “cắt và bỏ” không đánh giá được: cần khoảng theo dõi sau cắt polyp và mô bệnh học của polyp > 5 mm, không có trong bộ dữ liệu."))
ft312bTiêu chí | Bác sĩ nội soi | CADx | Đánh giá – bác sĩ | Đánh giá – CADx |
|---|---|---|---|---|
PIVI – để lại tại chỗ: NPV ≥ 90% | 91,3 (81,6–98,0) | 94,9 (88,1–100,0) | ≈ Đạt ước lượng điểm | ≈ Đạt ước lượng điểm |
SODA – để lại tại chỗ: Se ≥ 90% và Sp ≥ 80% | Se 88,6 (76,5–97,4) | Se 91,4 (81,6–100,0) | ✗ Chưa đạt (độ nhạy, độ đặc hiệu) | ✗ Chưa đạt (độ đặc hiệu) |
SODA – cắt và bỏ: Se ≥ 80% và Sp ≥ 80% | Se 94,6 (90,1–98,0) | Se 87,7 (82,0–93,0) | ✗ Chưa đạt (độ đặc hiệu) | ✗ Chưa đạt (độ đặc hiệu) |
PIVI – cắt và bỏ: đồng thuận khoảng theo dõi ≥ 90% | — | — | Không đánh giá được | Không đánh giá được |
% (KTC 95% bootstrap cụm theo bệnh nhân, B = 2000). Se: độ nhạy; Sp: độ đặc hiệu; NPV: giá trị tiên đoán âm – tính cho polyp tân sinh (u tuyến + SSL). | ||||
✓ Đạt: cận dưới KTC 95% ≥ ngưỡng; ≈ Đạt ước lượng điểm: ước lượng ≥ ngưỡng nhưng cận dưới KTC < ngưỡng; ✗ Chưa đạt: ước lượng < ngưỡng (ghi chỉ số chưa đạt). | ||||
Ngưỡng: ASGE PIVI (Rex và cs., Gastrointest Endosc 2011;73:419–22); ESGE SODA (Houwen và cs., Endoscopy 2022;54:88–99). Các ngưỡng quy định cho chẩn đoán với mức độ tin cậy cao; nghiên cứu không ghi nhận mức độ tin cậy nên đây là đối chiếu tham khảo. | ||||
PIVI “cắt và bỏ” không đánh giá được: cần khoảng theo dõi sau cắt polyp và mô bệnh học của polyp > 5 mm, không có trong bộ dữ liệu. | ||||
save_as_docx(ft312b, path = file.path(thu_muc_kq, "Bang_3_12b_Doi_chieu_nguong_PIVI_SODA.docx"))Biểu đồ khoảng tin cậy (dạng forest plot): mỗi chỉ số một hàng, chấm = ước lượng điểm, thanh ngang = KTC 95% bootstrap cụm; hai đường đứt đoạn đánh dấu ngưỡng SODA 80% và PIVI 90%; cột bên phải ghi p so sánh bác sĩ với CADx. Vẽ 2 hình: toàn bộ đại tràng (Biểu đồ 3.5) và trực tràng – sigma (Biểu đồ 3.6), mỗi hình có bản tiếng Việt và tiếng Anh.
MAU_NG <- c(BS = "#7B2D8E", CADx = "#4A90C8") # tím – bác sĩ; xanh – CADx
ve_nguong <- function(diem_x, ci_x, p_x, lang = "vi", n_lab = "") {
en <- lang == "en"
cs <- c("Se", "Sp", "PPV", "NPV", "Acc")
ten <- if (en) c(Se = "Sensitivity", Sp = "Specificity", PPV = "PPV", NPV = "NPV", Acc = "Accuracy")
else c(Se = "Độ nhạy", Sp = "Độ đặc hiệu", PPV = "PPV", NPV = "NPV", Acc = "Độ chính xác")
so <- if (en) function(x, k = 1) sprintf("%.*f", k, x) else so_vn
fp_ <- if (en) function(p) ifelse(p < 0.001, "p < 0.001", sprintf("p = %.3f", p))
else function(p) ifelse(p < 0.001, "p < 0,001", paste("p =", sub(".", ",", sprintf("%.3f", p), fixed = TRUE)))
k <- length(cs)
dd <- bind_rows(lapply(seq_along(cs), function(i) data.frame(
cs = cs[i], hang = k - i + 1, pp = c("BS", "CADx"),
est = c(diem_x[[paste0("BS_", cs[i])]], diem_x[[paste0("CADx_", cs[i])]]),
lo = ci_x[1, paste0(c("BS_", "CADx_"), cs[i])], hi = ci_x[2, paste0(c("BS_", "CADx_"), cs[i])]))) |>
mutate(y = hang + ifelse(pp == "BS", 0.2, -0.2), pp = factor(pp, levels = c("BS", "CADx")))
x_min <- max(0, floor(min(dd$lo) * 10) / 10 - 0.05)
nen <- data.frame(hang = seq(k, 1, by = -2)) # dải nền xám xen kẽ
hang_p <- data.frame(hang = k:1, nhan = fp_(p_x[cs]))
hang_t <- data.frame(hang = k:1, nhan = ten[cs])
nhan_pp <- if (en) c(BS = "Endoscopist", CADx = "CADx") else c(BS = "Bác sĩ nội soi", CADx = "CADx")
lr <- sprintf(if (en) "LR+ %s vs %s; LR− %s vs %s." else "LR+ %s vs %s; LR− %s vs %s.",
so(diem_x[["BS_LRp"]], 2), so(diem_x[["CADx_LRp"]], 2), so(diem_x[["BS_LRn"]], 2), so(diem_x[["CADx_LRn"]], 2))
chu <- if (en) sprintf("%s. 95%% CI: patient-level cluster bootstrap (B = %d).\np: GEE (sensitivity, specificity, accuracy) or cluster bootstrap (PPV, NPV). %s",
n_lab, B, lr)
else sprintf("%s. KTC 95%% bootstrap theo cụm bệnh nhân (B = %d).\np: GEE (độ nhạy, độ đặc hiệu, độ chính xác) hoặc bootstrap cụm (PPV, NPV). %s",
n_lab, B, lr)
ggplot(dd) +
geom_rect(data = nen, aes(xmin = -Inf, xmax = Inf, ymin = hang - 0.5, ymax = hang + 0.5),
fill = "#F3F4F6", inherit.aes = FALSE) +
geom_vline(xintercept = seq(0.3, 1, 0.1), colour = "grey88", linewidth = 0.3) +
geom_vline(xintercept = c(NGUONG_SODA, NGUONG_PIVI), linetype = "dashed", colour = "#1F4E5A", linewidth = 0.4) +
annotate("text", x = c(NGUONG_SODA, NGUONG_PIVI), y = k + 0.72, label = c("SODA 80%", "PIVI 90%"),
hjust = c(1.05, -0.05), size = 2.5, fontface = "bold", colour = "#1F4E5A", family = FONT) +
geom_segment(aes(x = lo, xend = hi, y = y, yend = y, colour = pp), linewidth = 0.9, show.legend = FALSE) +
geom_point(aes(x = est, y = y, fill = pp), shape = 21, colour = "white", size = 2.6, stroke = 0.6) +
geom_text(aes(x = hi, y = y, label = so(100 * est)), hjust = -0.25, size = 2.5, family = FONT, colour = "grey20") +
geom_text(data = hang_t, aes(x = x_min - 0.02, y = hang, label = nhan), hjust = 1, size = 2.9,
fontface = "bold", family = FONT, inherit.aes = FALSE) +
geom_text(data = hang_p, aes(x = 1.19, y = hang, label = nhan), hjust = 1, size = 2.7,
fontface = "bold", family = FONT, inherit.aes = FALSE) +
scale_colour_manual(values = MAU_NG, labels = nhan_pp) +
scale_fill_manual(values = MAU_NG, labels = nhan_pp,
guide = guide_legend(title = if (en) " ● point estimate — 95% CI" else " ● ước lượng điểm — KTC 95%",
title.position = "right", title.theme = element_text(size = 8, colour = "grey35", family = FONT))) +
scale_x_continuous(breaks = seq(0.3, 1, 0.1), labels = percent_format(accuracy = 1),
expand = expansion(mult = c(0, 0))) +
scale_y_continuous(breaks = NULL, expand = expansion(add = c(0.1, 0.4))) +
coord_cartesian(xlim = c(x_min, 1.19), clip = "off") +
labs(x = NULL, y = NULL, caption = chu,
subtitle = NULL) +
theme_bao() +
theme(axis.line = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(),
legend.justification = "left", plot.caption.position = "plot",
plot.margin = margin(6, 8, 4, 70))
}lab_all_vi <- sprintf("Toàn bộ đại tràng, %s", n_all)
lab_all_en <- sprintf("Whole colon, n = %d polyps/%d patients", nrow(d_pair), n_distinct(d_pair$id_bn))
g35_vi <- ve_nguong(diem, ci, p_chinh, "vi", lab_all_vi)
g35_en <- ve_nguong(diem, ci, p_chinh, "en", lab_all_en)
luu_hinh(g35_vi, "Bieu_do_3_5_Doi_chieu_nguong_toan_bo_VN", 17.5, 10.5)
luu_hinh(g35_en, "Figure_3_5_Thresholds_whole_colon_EN", 17.5, 10.5)
g35_vi
g35_en
lab_rs_vi <- sprintf("Trực tràng – sigma, %s", n_rs)
lab_rs_en <- sprintf("Rectosigmoid, n = %d polyps/%d patients", nrow(d_rs), n_distinct(d_rs$id_bn))
g36_vi <- ve_nguong(diem_rs, ci_rs, p_rs, "vi", lab_rs_vi)
g36_en <- ve_nguong(diem_rs, ci_rs, p_rs, "en", lab_rs_en)
luu_hinh(g36_vi, "Bieu_do_3_6_Doi_chieu_nguong_truc_trang_sigma_VN", 17.5, 10.5)
luu_hinh(g36_en, "Figure_3_6_Thresholds_rectosigmoid_EN", 17.5, 10.5)
g36_vi
g36_en
PHẦN 9. KẾT QUẢ KHÔNG XÁC ĐỊNH VÀ DỮ LIỆU THIẾU (STARD mục 15, 16, 25)
Báo cáo số polyp thiếu kết quả của từng phương pháp và phân tích độ nhạy theo kịch bản xấu nhất (coi mọi kết quả CADx thiếu là chẩn đoán sai). Biến cố bất lợi không có trong file dữ liệu nên phải điền thủ công từ hồ sơ.
d_wc <- d |> filter(!is.na(mbh_n), !is.na(bs_ts)) |>
mutate(cadx_wc = ifelse(is.na(cadx_ts), 1L - ref_ts, cadx_ts))
s_wc <- chi_so(d_wc$cadx_wc, d_wc$ref_ts)
tb313 <- data.frame(
`Nội dung` = c("Polyp CADx không cho kết quả, n", "Polyp thiếu chẩn đoán bác sĩ, n",
"Polyp thiếu mô bệnh học, n",
"Độ nhạy / độ đặc hiệu CADx – phân tích chính, %",
"Độ nhạy / độ đặc hiệu CADx – kịch bản xấu nhất (coi là sai), %",
"Biến cố bất lợi liên quan đến test chỉ điểm / tiêu chuẩn tham chiếu"),
`Kết quả` = c(sum(is.na(d$cadx_ts)), sum(is.na(d$bs_ts)), sum(is.na(d$mbh_n)),
sprintf("%.1f / %.1f", 100 * diem[["CADx_Se"]], 100 * diem[["CADx_Sp"]]),
sprintf("%.1f / %.1f", 100 * s_wc[["Se"]], 100 * s_wc[["Sp"]]),
"Ghi nhận từ hồ sơ (điền thủ công)"),
check.names = FALSE)
dinh_dang_bang(flextable(tb313) |> align(j = 2, align = "center", part = "all"),
"Bảng 3.13. Kết quả không xác định, dữ liệu thiếu và phân tích độ nhạy")Nội dung | Kết quả |
|---|---|
Polyp CADx không cho kết quả, n | 0 |
Polyp thiếu chẩn đoán bác sĩ, n | 0 |
Polyp thiếu mô bệnh học, n | 0 |
Độ nhạy / độ đặc hiệu CADx – phân tích chính, % | 87.7 / 62.1 |
Độ nhạy / độ đặc hiệu CADx – kịch bản xấu nhất (coi là sai), % | 87.7 / 62.1 |
Biến cố bất lợi liên quan đến test chỉ điểm / tiêu chuẩn tham chiếu | Ghi nhận từ hồ sơ (điền thủ công) |
PHẦN 10. XUẤT KẾT QUẢ
Ghi toàn bộ bảng vào một file Excel (mỗi bảng một sheet) và liệt kê các file đã tạo trong thư mục Output/ket_qua/.
writexl::write_xlsx(
list(Bang_0 = kiem_tra, Bang_3_1 = tb31, Bang_3_2 = tb32, Bang_3_3 = tb33,
Bang_3_4 = bang_2x2(d_pair$bs_ts, d_pair$ref_ts, "Bác sĩ"),
Bang_3_5 = bang_2x2(d_pair$cadx_ts, d_pair$ref_ts, "CADx"),
Bang_3_6 = tb36, Bang_3_6b = tbw,
Bang_3_7 = bang_2x2(d_pair$cadx_ts, d_pair$bs_ts, "CADx", "BS"), Bang_3_8 = tb38,
Bang_3_9 = bang3(d_pair$bs3, "Bác sĩ"), Bang_3_10 = bang3(d_pair$cadx3, "CADx"),
Bang_3_11 = tb311, Bang_3_12 = tb312,
Bang_3_12b = tb312b[, 1:5], Bang_3_13 = tb313),
file.path(thu_muc_kq, "Cac_bang_ket_qua_V2.xlsx"))
file_moi <- list.files(thu_muc_kq, pattern = "(_VN|_EN|V2|compareGroups|PIVI_SODA)\\.(png|tiff|xlsx|docx)$")
knitr::kable(data.frame(`File đã tạo trong Output/ket_qua/` = file_moi, check.names = FALSE))| File đã tạo trong Output/ket_qua/ |
|---|
| Bang_3_1_compareGroups.docx |
| Bang_3_1_compareGroups.xlsx |
| Bang_3_12b_Doi_chieu_nguong_PIVI_SODA.docx |
| Bang_3_2_compareGroups.docx |
| Bang_3_2_compareGroups.xlsx |
| Bieu_do_3_1_Vi_tri_VN.png |
| Bieu_do_3_1_Vi_tri_VN.tiff |
| Bieu_do_3_2_Phan_bo_chan_doan_VN.png |
| Bieu_do_3_2_Phan_bo_chan_doan_VN.tiff |
| Bieu_do_3_3_Chi_so_chan_doan_VN.png |
| Bieu_do_3_3_Chi_so_chan_doan_VN.tiff |
| Bieu_do_3_4_Ma_tran_nham_lan_VN.png |
| Bieu_do_3_4_Ma_tran_nham_lan_VN.tiff |
| Bieu_do_3_5_Doi_chieu_nguong_toan_bo_VN.png |
| Bieu_do_3_5_Doi_chieu_nguong_toan_bo_VN.tiff |
| Bieu_do_3_6_Doi_chieu_nguong_truc_trang_sigma_VN.png |
| Bieu_do_3_6_Doi_chieu_nguong_truc_trang_sigma_VN.tiff |
| Cac_bang_ket_qua_V2.xlsx |
| Figure_3_1_Location_EN.png |
| Figure_3_1_Location_EN.tiff |
| Figure_3_2_Diagnosis_distribution_EN.png |
| Figure_3_2_Diagnosis_distribution_EN.tiff |
| Figure_3_3_Diagnostic_performance_EN.png |
| Figure_3_3_Diagnostic_performance_EN.tiff |
| Figure_3_4_Confusion_matrix_EN.png |
| Figure_3_4_Confusion_matrix_EN.tiff |
| Figure_3_5_Thresholds_whole_colon_EN.png |
| Figure_3_5_Thresholds_whole_colon_EN.tiff |
| Figure_3_6_Thresholds_rectosigmoid_EN.png |
| Figure_3_6_Thresholds_rectosigmoid_EN.tiff |
sessionInfo() # ghi lại phiên bản R và các gói – phục vụ tính lặp lại (reproducibility)## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.6.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Asia/Ho_Chi_Minh
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] scales_1.4.0 patchwork_1.3.2 ggsignif_0.6.4
## [4] ggpubr_1.0.0 ggplot2_4.0.3 DTComPair_1.2.6
## [7] PropCIs_0.3-0 geepack_1.3.13 officer_0.7.6
## [10] flextable_0.10.1 table1_1.5.1 compareGroups_4.10.4
## [13] tidyr_1.3.2 dplyr_1.2.1 readxl_1.5.0.1
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.6 gee_4.13-30 writexl_2.0.1
## [4] rlang_1.3.0 magrittr_2.0.5 PMCMRplus_1.9.12
## [7] otel_0.2.0 compiler_4.6.1 BWStest_0.2.3
## [10] systemfonts_1.3.2 vctrs_0.7.3 kSamples_1.2-12
## [13] stringr_1.6.0 pkgconfig_2.0.3 shape_1.4.6.1
## [16] fastmap_1.2.0 backports_1.5.1 labeling_0.4.3
## [19] rmdformats_1.0.4 utf8_1.2.6 rmarkdown_2.32
## [22] nloptr_2.2.1 ragg_1.5.2 purrr_1.2.2
## [25] xfun_0.61 glmnet_5.1 jomo_2.7-6
## [28] cachem_1.1.0 Rmpfr_1.1-3 jsonlite_2.0.0
## [31] SuppDists_1.1-9.9 gmp_0.7-5.1 uuid_1.2-2
## [34] pan_2.0 broom_1.0.13 parallel_4.6.1
## [37] R6_2.6.1 bslib_0.12.0 stringi_1.8.9
## [40] RColorBrewer_1.1-3 Rsolnp_2.0.1 car_3.1-5
## [43] parallelly_1.48.0 boot_1.3-32 rpart_4.1.27
## [46] jquerylib_0.1.4 cellranger_1.1.0 numDeriv_2016.8-1.1
## [49] assertthat_0.2.1 bookdown_0.48 Rcpp_1.1.2
## [52] iterators_1.0.14 knitr_1.52 future.apply_1.20.2
## [55] Matrix_1.7-5 splines_4.6.1 nnet_7.3-20
## [58] tidyselect_1.2.1 abind_1.4-8 rstudioapi_0.19.0
## [61] yaml_2.3.12 codetools_0.2-20 listenv_1.0.0
## [64] lattice_0.22-9 tibble_3.3.1 withr_3.0.3
## [67] S7_0.2.2 askpass_1.2.1 evaluate_1.0.5
## [70] future_1.76.0 survival_3.8-6 zip_3.0.2
## [73] xml2_1.6.0 pillar_1.11.1 carData_3.0-6
## [76] mice_3.19.0 foreach_1.5.2 reformulas_0.4.4
## [79] generics_0.1.4 truncnorm_1.0-9 minqa_1.2.8
## [82] chron_2.3-63 globals_0.19.1 glue_1.8.1
## [85] gdtools_0.5.1 tools_4.6.1 data.table_1.18.6.1
## [88] lme4_2.0-6 mvtnorm_1.4-2 grid_4.6.1
## [91] rbibutils_2.4.1 nlme_3.1-169 Formula_1.2-6
## [94] cli_3.6.6 kableExtra_1.4.1 textshaping_1.0.5
## [97] fontBitstreamVera_0.1.1 viridisLite_0.4.3 svglite_2.2.2
## [100] gtable_0.3.6 rstatix_1.1.0 sass_0.4.10
## [103] digest_0.6.39 fontquiver_0.2.1 farver_2.1.2
## [106] memoise_2.0.1 htmltools_0.5.9 lifecycle_1.0.5
## [109] multcompView_0.1-12 mitml_0.4-5 fontLiberation_0.1.0
## [112] openssl_2.4.2 MASS_7.3-65 HardyWeinberg_1.7.9 ```