Trong R, hồi quy tuyến tính chạy bằng lm(), hồi quy logistic bằng glm(..., family = binomial). Sau khi kiểm giả định, dựng bảng bằng gtsummary::tbl_regression() hoặc modelsummary(), chuyển sang flextable rồi save_as_docx() là có tệp Word, không phải chép số bằng tay.
Chép số từ cửa sổ kết quả sang Word là chỗ sinh lỗi nhiều nhất trong luận văn: đảo dấu, nhầm dòng, quên cập nhật khi chạy lại mô hình với dữ liệu đã sửa. Đưa bảng ra thẳng từ mã giúp mỗi lần chạy lại là bảng tự đúng theo. Ví dụ dưới đây dùng bộ dữ liệu có sẵn trong R để bạn chạy thử ngay; mã đã chạy với gtsummary 2.6, modelsummary 2.6 và flextable 0.10.
Chạy mô hình và đọc kết quả gốc
Mô hình tuyến tính: m1 <- lm(mpg ~ wt, data = mtcars) rồi m2 <- lm(mpg ~ wt + hp + am, data = mtcars). Lệnh summary(m2) in hệ số, sai số chuẩn, giá trị t, p, R bình phương và R bình phương hiệu chỉnh. Khoảng tin cậy lấy bằng confint(m2).
Biến phân loại nên đổi sang factor trước, và chọn nhóm tham chiếu có chủ đích bằng relevel(x, ref = 'Nhóm chuẩn'). Nhóm tham chiếu mặc định là mức đầu tiên theo thứ tự chữ cái, hiếm khi là nhóm bạn muốn so.
Mô hình logistic: m3 <- glm(vs ~ wt + hp, data = mtcars, family = binomial). Hệ số in ra ở thang log odds; lấy tỷ số chênh bằng exp(coef(m3)). Trong bảng báo cáo, thường trình bày tỷ số chênh kèm khoảng tin cậy.
Hai việc nên làm ngay sau khi chạy. Thứ nhất, xem nobs(m2): lm() và glm() mặc định bỏ mọi dòng thiếu ở bất kỳ biến nào trong công thức, nên số quan sát thật dùng trong mô hình có thể nhỏ hơn cỡ mẫu nhiều. Khi so hai mô hình lồng nhau, cả hai phải chạy trên cùng một tập dòng; lọc trước bằng tidyr::drop_na() trên đủ các biến của mô hình lớn. Thứ hai, so mô hình bằng anova(m1, m2), cho kiểm định F xem các biến thêm vào có cải thiện mô hình không; với glm(), dùng anova(m_nho, m_lon, test = 'Chisq').
Kiểm giả định nào, bằng lệnh gì?
Giả định không đạt chưa chắc là thảm hoạ, nhưng phải biết để chọn cách xử lý và ghi vào bài:
| Giả định | Cách kiểm trong R | Nếu không đạt |
|---|---|---|
| Quan hệ tuyến tính | Hình phần dư theo giá trị dự báo: plot(m2, which = 1) | Thêm biến bậc hai, biến đổi biến, hoặc dùng spline |
| Phương sai sai số không đổi | plot(m2, which = 3), performance::check_heteroscedasticity(m2) | Dùng sai số chuẩn bền HC3 |
| Phần dư gần phân phối chuẩn | Hình Q-Q: plot(m2, which = 2) | Mẫu lớn thường ít đáng lo; mẫu nhỏ cân nhắc bootstrap |
| Đa cộng tuyến | car::vif(m2) | Xem lại biến trùng ý nghĩa, gộp hoặc bỏ bớt |
| Điểm ảnh hưởng lớn | plot(m2, which = 5), khoảng cách Cook | Kiểm lỗi nhập liệu, chạy thêm mô hình không có điểm đó để so |
| Tổng hợp nhanh | performance::check_model(m2) | Đọc từng khung, đừng chỉ nhìn chữ OK |
Với sai số chuẩn bền, cách gọn là truyền thẳng vào bảng: modelsummary(m2, vcov = 'HC3'). Nếu muốn xem riêng: lmtest::coeftest(m2, vcov = sandwich::vcovHC(m2, type = 'HC3')). Ghi rõ trong chú thích bảng là đã dùng sai số chuẩn bền loại nào.
Dựng bảng bằng gtsummary
Một mô hình, nhãn tiếng Việt
bang <- gtsummary::tbl_regression(m2, label = list(wt ~ 'Khối lượng xe', hp ~ 'Công suất')). Hàm tự đổi tên cột, gộp khoảng tin cậy thành một cột và đánh dấu nhóm tham chiếu. Thêm add_glance_table(include = c(r.squared, nobs)) để có R bình phương và số quan sát ở cuối bảng; hàm này cần gói broom.helpers, R sẽ nhắc nếu thiếu.
Logistic với tỷ số chênh
gtsummary::tbl_regression(m3, exponentiate = TRUE) in tỷ số chênh thay cho hệ số log odds, khoảng tin cậy cũng được lấy mũ theo.
Ghi ra Word
bang |> gtsummary::as_flex_table() |> flextable::save_as_docx(path = here::here('ket-qua', 'bang', 'bang-2.docx')). Mở tệp Word ra là thấy bảng, chép vào bản thảo hoặc chèn bằng cách dán.
Nhiều mô hình cạnh nhau bằng modelsummary
Bảng so sánh mô hình lồng nhau, mỗi mô hình một cột, là kiểu bảng quen trong bài báo kinh tế và xã hội học. modelsummary(list('Mô hình 1' = m1, 'Mô hình 2' = m2), stars = TRUE, gof_map = c('nobs', 'r.squared', 'adj.r.squared'), output = 'flextable') trả về một flextable, nối tiếp |> flextable::save_as_docx(path = ...).
Có thể ghi thẳng output = 'bang.docx', nhưng cách đó cần thêm công cụ pandoc trên máy; khi chạy thử, R báo yêu cầu cài thêm. Đi qua flextable tránh được phụ thuộc này.
Đổi tên biến cho cả bảng bằng coef_rename = c(wt = 'Khối lượng xe', hp = 'Công suất'). Muốn bỏ hệ số chặn khỏi bảng, dùng coef_omit = 'Intercept'.
Diễn giải hệ số cho đúng
Hệ số chưa chuẩn hoá và đơn vị
Trong m2, hệ số của wt cho biết mpg thay đổi bao nhiêu khi khối lượng xe tăng một đơn vị (ở bộ này là một nghìn pound), giữ các biến khác cố định. Diễn giải luôn đi kèm đơn vị của cả hai biến; nếu đơn vị gốc khó hình dung, đổi đơn vị trước khi chạy mô hình, ví dụ thu nhập tính theo triệu đồng thay vì đồng.
Hệ số chuẩn hoá
Một số ngành quen báo cáo beta chuẩn hoá. Cách minh bạch là chuẩn hoá biến trước bằng scale() rồi chạy lại mô hình, hoặc dùng effectsize::standardize_parameters(m2). Ghi rõ đã chuẩn hoá biến nào; với biến nhị phân, chuẩn hoá thường làm hệ số khó hiểu hơn.
Biến tương tác
Viết lm(mpg ~ wt * am, data = mtcars) là có cả hai tác động chính và tích wt:am. Khi có tương tác, hệ số của wt chỉ là tác động của khối lượng ở nhóm tham chiếu của am, không phải tác động chung. Trung tâm hoá biến liên tục (trừ trung bình) giúp hệ số chính dễ đọc hơn; vẽ đường dự báo theo nhóm, ví dụ bằng gói marginaleffects hoặc ggeffects, thường giúp người đọc hiểu nhanh hơn một bảng hệ số.
Logistic: tỷ số chênh không phải tỷ số nguy cơ
Tỷ số chênh 1,5 không có nghĩa xác suất tăng 50%. Khi kết cục phổ biến, khác biệt giữa hai khái niệm đáng kể. Muốn nói bằng xác suất, tính xác suất dự báo hoặc tác động biên trung bình và báo cáo kèm.
Khớp bảng với mẫu Word của trường hoặc tạp chí
Nhiều trường yêu cầu phông Times New Roman cỡ 13 cho thân luận văn và cỡ nhỏ hơn cho bảng. Đặt một lần ở đầu kịch bản để mọi bảng theo cùng một kiểu:
flextable::set_flextable_defaults(font.family = 'Times New Roman', font.size = 11, padding = 3)áp cho mọi bảng tạo sau đó.- Muốn bảng nằm trong một tệp có sẵn kiểu chữ và lề của trường:
officer::read_docx(here::here('tai-lieu', 'mau.docx')) |> flextable::body_add_flextable(ft) |> print(target = 'ket-qua.docx'). - Chú thích dưới bảng (nguồn, kiểu sai số chuẩn, mức ý nghĩa) thêm bằng
flextable::add_footer_lines(ft, 'Ghi chú: ...'). - Số thập phân: thống nhất hai hoặc ba chữ số trong toàn bài; dấu thập phân là dấu chấm hay phẩy theo yêu cầu nơi nộp.
Khi nộp tạp chí quốc tế, theo mẫu bảng trong hướng dẫn tác giả của chính tạp chí đó; khi nộp luận văn, theo quy định trình bày của cơ sở đào tạo.
Lỗi hay gặp khi báo cáo hồi quy
- Báo cáo p mà thiếu hệ số và khoảng tin cậy. Người đọc cần biết độ lớn của tác động, không chỉ biết có ý nghĩa hay không.
- Để hệ số log odds trong bảng logistic nhưng diễn giải như tỷ số chênh.
- Không ghi số quan sát thực dùng trong mô hình;
lm()bỏ dòng thiếu ở bất kỳ biến nào, nên số này thường nhỏ hơn cỡ mẫu khảo sát. - Dán bảng xong rồi sửa số tay trong Word khi phát hiện lỗi dữ liệu. Sửa trong mã và xuất lại.
- Thêm biến kiểm soát vì “ai cũng thêm” mà không giải thích vai trò của chúng; người phản biện sẽ hỏi vì sao có biến này, không có biến kia.
- Đọc R bình phương thấp như dấu hiệu mô hình sai. Trong khoa học xã hội, mô hình giải thích được một phần nhỏ phương sai vẫn có thể cho ước lượng hữu ích về một quan hệ cụ thể.
Câu hỏi thường gặp
Nên chọn gtsummary hay modelsummary?
gtsummary mạnh ở bảng một mô hình kiểu y sinh, có sẵn tỷ số chênh và bảng mô tả mẫu. modelsummary tiện cho bảng nhiều mô hình cạnh nhau và đổi kiểu sai số chuẩn. Có thể dùng cả hai trong cùng dự án.
Bảng xuất ra có giữ được tiếng Việt không?
Có. flextable ghi tệp Word theo Unicode nên nhãn tiếng Việt hiện đúng, miễn là phông chữ chọn có đủ dấu tiếng Việt, như Times New Roman hay Arial.
Làm sao xuất bảng thống kê mô tả mẫu ra Word?
Dùng gtsummary::tbl_summary trên bảng dữ liệu, chọn biến cần mô tả, chia nhóm bằng tham số by, rồi as_flex_table và save_as_docx như bảng hồi quy.
Có cần báo cáo kết quả kiểm giả định trong bài không?
Thường chỉ cần một hai câu trong phần phương pháp nói đã kiểm những gì và xử lý ra sao. Hình chẩn đoán có thể đưa vào phụ lục hoặc tài liệu bổ sung nếu tạp chí cho phép.
Biến phụ thuộc là thang Likert một câu thì dùng lm được không?
Một câu Likert là biến thứ bậc, mô hình logistic thứ bậc như MASS::polr hoặc gói ordinal thường hợp hơn. Điểm trung bình của cả thang đo nhiều câu thì hay được phân tích bằng lm.
