Hồi quy tuyến tính trong Stata là kỹ năng nền tảng của hầu hết luận văn kinh tế, quản trị và y tế công cộng. Chỉ với vài dòng do-file Stata, bạn ước lượng được mô hình, kiểm tra đa cộng tuyến bằng lệnh vif, kiểm định phương sai sai số thay đổi bằng estat hettest và chuyển sang sai số chuẩn bền vững (robust) khi cần. Bài viết này đi từng bước, giải thích từng con số trong cửa sổ Results và gợi ý cách trình bày kết quả trong luận văn. Ví dụ dùng một bộ dữ liệu mô phỏng 400 quan sát về sự hài lòng của khách hàng, chạy bằng công cụ Stata trực tuyến của ResearchSTAT.
Chuẩn bị dữ liệu và câu hỏi nghiên cứu
Giả sử đề tài muốn kiểm định ba giả thuyết: chất lượng dịch vụ (CLD), giá cả cảm nhận (GCA) và niềm tin (NIT) tác động cùng chiều đến sự hài lòng (SHL). Mỗi khái niệm đo bằng 4 câu hỏi Likert 5 điểm. Trước khi chạy hồi quy, bạn tạo biến đại diện là trung bình các câu hỏi của từng thang đo, sau khi đã kiểm tra độ tin cậy bằng Cronbach's Alpha (xem bài Cronbach's Alpha). Trong Stata, biến trung bình được tạo bằng lệnh generate:
* Tạo biến đại diện cho từng thang đo
generate CLD = (CLD1 + CLD2 + CLD3 + CLD4) / 4
generate SHL = (SHL1 + SHL2 + SHL3 + SHL4) / 4
describe CLD SHL
summarize CLD GCA NIT SHLLưu ý: trong Stata, phép cộng có một giá trị thiếu sẽ cho kết quả thiếu. Nếu muốn lấy trung bình của các câu còn lại, dùng egen CLD = rowmean(CLD1-CLD4). Lệnh summarize cho số quan sát (Obs), trung bình, độ lệch chuẩn, giá trị nhỏ nhất và lớn nhất; đây là bước để phát hiện giá trị nhập sai ngoài khoảng 1-5 trước khi phân tích.
Chạy lệnh regress và đọc phần đầu bảng kết quả
Cú pháp cơ bản là regress biến_phụ_thuộc danh_sách_biến_độc_lập. Với ví dụ của chúng ta:
regress SHL CLD GCA NITPhần đầu bảng gồm hai khối. Khối bên trái là bảng phân tích phương sai: tổng bình phương của mô hình (Model SS = 79,586), của phần dư (Residual SS = 173,976) và tổng (Total SS = 253,562). Khối bên phải là các chỉ số tóm tắt:
- Number of obs = 400: số quan sát thực sự được dùng. Nếu nhỏ hơn cỡ mẫu, đã có quan sát bị loại vì thiếu dữ liệu.
- F(3, 396) = 60,38 và Prob > F = 0,0000: kiểm định F cho giả thuyết mọi hệ số hồi quy đồng thời bằng 0. Prob > F nhỏ hơn 0,05 nghĩa là mô hình có ý nghĩa thống kê.
- R-squared = 0,3139: ba biến độc lập giải thích 31,39% biến thiên của sự hài lòng.
- Adj R-squared = 0,3087: R² hiệu chỉnh theo số biến, dùng khi so sánh các mô hình có số biến khác nhau.
- Root MSE = 0,663: độ lệch chuẩn của phần dư, đơn vị là đơn vị của biến phụ thuộc.
Đọc bảng hệ số: Coefficient, Std. err., t, P>|t|
Mỗi dòng của bảng hệ số tương ứng với một biến độc lập, dòng cuối _cons là hệ số chặn.
| Biến | Coefficient | Std. err. | t | Giá trị p | KTC 95% |
|---|---|---|---|---|---|
| CLD | 0,238 | 0,043 | 5,58 | 0,000 | [0,154; 0,322] |
| GCA | 0,183 | 0,043 | 4,24 | 0,000 | [0,098; 0,267] |
| NIT | 0,393 | 0,046 | 8,52 | 0,000 | [0,303; 0,484] |
| _cons | 0,608 | 0,214 | 2,85 | 0,005 | [0,188; 1,028] |
Cách đọc: giữ các biến khác không đổi, khi niềm tin tăng 1 điểm thì sự hài lòng trung bình tăng 0,393 điểm. Giá trị t là hệ số chia cho sai số chuẩn; P>|t| là giá trị p của kiểm định hai phía. Cả ba biến đều có p < 0,001 nên ba giả thuyết được ủng hộ. Khoảng tin cậy 95% không chứa 0 cũng cho cùng kết luận. Nếu muốn so sánh mức độ tác động giữa các biến có đơn vị khác nhau, thêm tùy chọn beta để Stata in hệ số chuẩn hóa:
regress SHL CLD GCA NIT, betaKiểm tra đa cộng tuyến với lệnh vif
Đa cộng tuyến xảy ra khi các biến độc lập tương quan mạnh với nhau, làm sai số chuẩn phình to và hệ số đổi dấu bất thường. Sau regress, gõ vif (hoặc estat vif) để xem hệ số phóng đại phương sai của từng biến.
Trong ví dụ, VIF của CLD và NIT bằng 1,16, của GCA bằng 1,00, trung bình 1,11. Nguyên tắc thực hành phổ biến: VIF dưới 5 là chấp nhận được, trên 10 là đa cộng tuyến nghiêm trọng (Hair và cộng sự, 2019). Với thang đo Likert, VIF cao thường do hai khái niệm đo trùng nội dung; khi đó hãy xem lại kết quả EFA hoặc kiểm tra giá trị phân biệt thay vì loại biến một cách máy móc.
Kiểm định phương sai sai số thay đổi với estat hettest
Hồi quy bình phương nhỏ nhất (OLS) giả định phương sai của sai số không đổi. Khi giả định bị vi phạm, hệ số vẫn không chệch nhưng sai số chuẩn sai, kéo theo giá trị p không còn đáng tin. Hai cách kiểm tra bổ sung cho nhau: vẽ đồ thị phần dư theo giá trị dự báo bằng rvfplot và chạy kiểm định Breusch-Pagan/Cook-Weisberg bằng estat hettest (Breusch và Pagan, 1979; Cook và Weisberg, 1983).
Kết quả estat hettest trong ví dụ: chi2(1) = 1,81 và Prob > chi2 = 0,1789. Vì giá trị p lớn hơn 0,05, chưa đủ cơ sở bác bỏ giả thuyết H0 phương sai không đổi. Nếu p nhỏ hơn 0,05, đồ thị phần dư thường có dạng phễu: sai số lớn dần theo giá trị dự báo.
Khi giả định bị vi phạm: sai số chuẩn bền vững
Cách xử lý phổ biến nhất là ước lượng lại với sai số chuẩn bền vững theo White, trong Stata là tùy chọn robust hoặc vce(robust) (dạng HC1). Hệ số không thay đổi, chỉ sai số chuẩn, t, p và khoảng tin cậy được tính lại.
regress SHL CLD GCA NIT, robust
Trong ví dụ, sai số chuẩn robust của CLD là 0,046 so với 0,043 khi chưa hiệu chỉnh; kết luận về ý nghĩa thống kê không đổi. Nhiều tạp chí kinh tế khuyến khích báo cáo sai số chuẩn robust ngay từ đầu. Nếu dữ liệu có cấu trúc nhóm (nhiều nhân viên trong cùng một chi nhánh), cân nhắc vce(cluster tên_biến_nhóm).
Quy trình hồi quy trong một do-file
Một do-file tốt giúp kết quả tái lập được: người hướng dẫn hoặc phản biện có thể chạy lại và nhận đúng các con số trong luận văn. Gợi ý trình tự như sau.
* 1. Mô tả dữ liệu
describe
summarize CLD GCA NIT SHL
* 2. Tương quan
pwcorr CLD GCA NIT SHL, sig obs
* 3. Hồi quy
regress SHL CLD GCA NIT
* 4. Chẩn đoán
vif
estat hettest
* 5. Nếu cần: sai số chuẩn bền vững
regress SHL CLD GCA NIT, robustTrình bày kết quả hồi quy trong luận văn
Một đoạn trình bày mẫu theo kiểu APA:
Kết quả hồi quy tuyến tính bội cho thấy mô hình có ý nghĩa thống kê (F(3, 396) = 60,38; p < 0,001) và giải thích 31,4% biến thiên của sự hài lòng (R² = 0,314; R² hiệu chỉnh = 0,309). Niềm tin (B = 0,393; p < 0,001), chất lượng dịch vụ (B = 0,238; p < 0,001) và giá cả cảm nhận (B = 0,183; p < 0,001) đều tác động cùng chiều đến sự hài lòng. Không có dấu hiệu đa cộng tuyến (VIF ≤ 1,16) và phương sai sai số thay đổi (Breusch-Pagan χ²(1) = 1,81; p = 0,179).
Ghi rõ phần mềm đã dùng trong mục phương pháp. Nếu bạn chạy trên ResearchSTAT, hãy ghi "phân tích bằng trình thông dịch cú pháp Stata của ResearchSTAT" thay vì ghi Stata, vì đây là bộ máy tính toán độc lập, không phải phần mềm Stata của StataCorp. Kết quả trên cùng dữ liệu đã được đối chiếu khớp với hàm lm() của R (xem bài phân tích dữ liệu bằng R).
Chạy thử ngay trên ResearchSTAT
Công cụ Stata của ResearchSTAT (giới thiệu chi tiết ở trang chạy lệnh Stata trực tuyến) cho phép bạn viết do-file, chạy với dữ liệu đã tải lên, xem kết quả theo bố cục cửa sổ Results, vẽ histogram, scatter, graph box, sau đó tải file .log, Word hoặc ảnh PNG. Lệnh được hỗ trợ gồm describe, summarize, tabulate, correlate, pwcorr, ttest, oneway, regress, logit, probit, alpha, vif, estat hettest, predict, generate, replace. Xem thêm bài hồi quy logistic và probit nếu biến phụ thuộc chỉ có hai giá trị.
Lệnh regress trong Stata có tự loại quan sát thiếu không?
Có. Stata dùng cách loại theo danh sách (listwise): mọi quan sát thiếu ở bất kỳ biến nào trong mô hình đều bị loại. Kiểm tra Number of obs để biết số quan sát thực dùng.
Nên báo cáo hệ số chuẩn hóa hay chưa chuẩn hóa?
Báo cáo hệ số chưa chuẩn hóa (Coefficient) kèm sai số chuẩn và p; thêm hệ số chuẩn hóa (tùy chọn beta) khi cần so sánh mức tác động giữa các biến.
VIF bao nhiêu thì phải xử lý?
Thường VIF trên 5 cần xem xét và trên 10 là nghiêm trọng. Với dữ liệu khảo sát, nên kiểm tra lại cách đo các khái niệm trước khi loại biến.
Kết quả trên ResearchSTAT có phải do phần mềm Stata tạo ra không?
Không. ResearchSTAT có trình thông dịch cú pháp Stata riêng, tính bằng numpy, scipy và statsmodels, trình bày theo bố cục Stata. Hãy ghi đúng công cụ trong luận văn.
