Hồi quy logistic trong R: glm() với family = binomial

Khi kết quả là có/không

Hồi quy tuyến tính dự đoán một con số. Nhưng nhiều câu hỏi đáng mô hình hóa lại là nhị phân: khách hàng có rời bỏ không, bệnh nhân có hồi phục không, email có được nhấp không. Khớp một đường thẳng với kết quả 0/1 hỏng ngay lập tức: đường thẳng vô tư dự đoán xác suất −0,3 hoặc 1,4, những con số vô nghĩa.

Hồi quy logistic sửa điều này bằng cách mô hình hóa xác suất của kết quả thông qua phép biến đổi log-odds (logit): log(p / (1 − p)) = intercept + slope × x. Thang log-odds trải dài cả trục số, nên một phương trình tuyến tính khớp tự nhiên ở đó, và khi chuyển ngược lại, mọi dự đoán bị ép vào khoảng (0, 1) theo đường cong hình chữ S quen thuộc. Cái giá của mẹo này: các hệ số nằm trên thang log-odds, và toàn bộ việc đọc một mô hình hồi quy logistic là dịch chúng ngược về thứ con người hiểu được.

Khớp mô hình: glm() với family = binomial

glm() (generalized linear model, mô hình tuyến tính tổng quát) là người anh lớn của lm(); family = binomial chọn hồi quy logistic. Trong mtcars, am ghi loại hộp số (1 = số sàn, 0 = số tự động): xe tiết kiệm nhiên liệu có xu hướng là số sàn không?

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(fit)

Hai điều trước khi đọc đầu ra. Thứ nhất, family = binomial không phải tùy chọn: bỏ nó đi, glm() sẽ lặng lẽ khớp bình phương nhỏ nhất thông thường. Thứ hai, biến kết quả phải là nhị phân: 0/1, logic, hoặc một factor hai mức (R mô hình hóa xác suất của mức thứ hai).

Giờ là summary, từng khối một:

• Coefficients: Estimate của mpg khoảng 0,31, và đó là độ dốc log-odds: mỗi mpg tăng thêm cộng 0,31 vào log-odds của việc là xe số sàn. Dương nghĩa là "làm tăng xác suất", âm nghĩa là "làm giảm"; ngoài dấu ra, trực giác của không ai hoạt động trên thang này, đó là lý do có phần tiếp theo.
• z value và Pr(>|z|): cùng logic với kiểm định t của hồi quy (Estimate ÷ Std. Error, rồi một p-value cho câu hỏi "có thể bằng 0 không?"), chỉ là dùng xấp xỉ chuẩn, nên là z thay vì t. Ở đây p ≈ 0,011: mối liên hệ giữa mpg và loại hộp số khó mà chỉ là nhiễu.
• Null deviance và Residual deviance: deviance là thước đo độ khớp kém trong thế giới glm (càng nhỏ càng tốt). Null deviance (43,2 với 31 bậc tự do) là của mô hình chỉ có hệ số chặn; residual deviance (29,7 với 30 bậc tự do) là của mô hình bạn. Mức giảm khoảng 13,6 đổi lấy một bậc tự do của biến dự báo là phiên bản glm của "R bình phương đã tăng".
• AIC: điểm so sánh mô hình, cân bằng độ khớp với độ phức tạp; thấp hơn thì thắng. Đứng một mình thì vô nghĩa, hữu ích khi so giữa các mô hình ứng viên trên cùng dữ liệu.

Từ log-odds sang tỷ số odds: exp(coef())

Lấy mũ chuyển các hệ số từ thang log-odds dạng cộng sang thang odds dạng nhân:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

print(exp(coef(fit)))
print(exp(confint(fit)))

exp(0.307) ≈ 1.36, và đây là mẫu câu trung thực nên thuộc lòng: "mỗi mpg tăng thêm nhân odds của xe số sàn lên khoảng 1,36 lần." Tỷ số odds lớn hơn 1 làm tăng odds, nhỏ hơn 1 làm giảm, đúng bằng 1 là không có tác động, đó là lý do kết luận từ khoảng tin cậy cho tỷ số odds là "khoảng có loại trừ 1 không?" (không phải 0; 0 là ranh giới khi còn ở thang log-odds).

Cẩn thận với cách diễn đạt: odds không phải xác suất. Odds = p / (1 − p), nên xác suất 0,75 là odds bằng 3. Nhân odds với 1,36 không giống nhân xác suất với 1,36, và khi kết quả phổ biến, khoảng cách đó rất lớn. Tỷ số odds bằng 2 cho một kết quả hiếm hoạt động gần giống "nguy cơ tăng khoảng gấp đôi"; cho một kết quả xảy ra ở mức 50% thì hoàn toàn không phải vậy. Đừng bao giờ báo cáo tỷ số odds bằng cách nói của tỷ số nguy cơ ("khả năng cao gấp 1,36 lần") trừ khi kết quả hiếm.

Xác suất dự đoán: cái bẫy type = "response"

Bug hồi quy logistic phổ biến nhất ngoài đời thực:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
new_cars <- data.frame(mpg = c(15, 20, 25, 30))

print(predict(fit, newdata = new_cars))                      # log-odds!
print(predict(fit, newdata = new_cars, type = "response"))   # probabilities

Lệnh gọi thứ nhất trả về theo mặc định type = "link": dự đoán trên thang log-odds, có cả giá trị âm. Lệnh gọi thứ hai trả về xác suất thật. Nếu "xác suất" của bạn có lúc âm hoặc lớn hơn 1, đây chính là lý do. Chạy khối code: một xe 15 mpg gần như không có khả năng là số sàn, một xe 30 mpg rất có thể là số sàn, và đường cong chữ S uốn qua khoảng giữa.

Phân loại: đặt ngưỡng và bảng nhầm lẫn

Xác suất trở thành lớp dự đoán khi bạn chọn một ngưỡng cắt (0,5 là lựa chọn mặc định), và bảng điểm trung thực là bảng dự đoán so với thực tế:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

probs <- predict(fit, type = "response")   # fitted probabilities, one per car
predicted <- probs > 0.5

print(table(actual = mtcars$am, predicted = predicted))

Các ô trên đường chéo là dự đoán đúng; hai ô ngoài đường chéo là hai loại sai lầm khác nhau (dự đoán số sàn cho xe số tự động, và ngược lại). Chỉ riêng độ chính xác tổng thể có thể tô vẽ quá mức cho một mô hình: nếu 95% khách hàng không rời bỏ, cách "dự đoán không ai rời bỏ" đạt 95% trong khi không bắt được ai rời bỏ cả, nên hãy luôn xem cả hai loại lỗi. Và 0,5 là quy ước, không phải luật: khi hai loại sai lầm có chi phí khác nhau, hãy dời ngưỡng cho phù hợp.

Một lưu ý trung thực: bảng này chấm điểm mô hình trên chính dữ liệu đã dùng để khớp, nên kết quả đẹp hơn thực tế. Đánh giá thật sự cần giữ lại dữ liệu mà mô hình chưa từng thấy.

Nhiều biến dự báo

Y như lm(): thêm số hạng bằng +, và mọi cách diễn giải được thêm điều kiện "khi giữ nguyên các biến còn lại":

fit2 <- glm(am ~ mpg + qsec, data = mtcars, family = binomial)
summary(fit2)
print(exp(coef(fit2)))

Giờ mỗi hệ số đã lấy mũ là hệ số nhân odds khi biến dự báo đó tăng một đơn vị giữa những xe giống nhau ở các biến dự báo khác. Cơ chế mở rộng được, nhưng các lưu ý từ hồi quy tuyến tính cũng vậy: các biến dự báo tương quan với nhau sẽ xáo trộn hệ số của nhau.

Những điều cần cẩn thận

• Phân tách hoàn toàn. Nếu một biến dự báo chia biến kết quả một cách hoàn hảo (mọi xe trên một mức mpg nào đó đều là số sàn, mọi xe dưới mức đó đều là số tự động), hệ số hợp lý cực đại sẽ muốn tiến ra vô cùng. R cảnh báo glm.fit: fitted probabilities numerically 0 or 1 occurred và báo cáo các hệ số khổng lồ với sai số chuẩn vô lý. Đừng dùng những con số đó; hãy đơn giản hóa mô hình, thu thêm dữ liệu, hoặc dùng phương pháp có phạt (các package brglm2 hoặc logistf).
• Đủ số sự kiện. Ràng buộc thật sự là số lượng của kết quả hiếm hơn, không phải tổng số hàng. Một quy tắc kinh nghiệm cũ cần khoảng từ 10 đến 15 sự kiện cho mỗi biến dự báo; các ví dụ 32 xe ở đây để dạy cơ chế, không phải khuôn mẫu cho cỡ mẫu đủ để công bố.
• Tỷ số odds không phải tỷ số nguy cơ khi kết quả phổ biến: đã nói ở trên, nhắc lại vì người phản biện sẽ phát hiện ra dù bạn không nhận thấy.

Những điều cần nhớ

• Kết quả nhị phân → glm(y ~ x, data = df, family = binomial); đừng bao giờ quên family.
• Hệ số thô là log-odds; exp(coef(fit)) cho tỷ số odds, và giá trị không có tác động cho khoảng tin cậy của chúng là 1.
• Mẫu câu: "mỗi lần x tăng một đơn vị nhân odds của kết quả lên exp(b) lần."
• predict(..., type = "response") để lấy xác suất: mặc định trả về log-odds, nhầm lẫn số một.
• Phân loại bằng một ngưỡng và đánh giá bằng bảng nhầm lẫn; chỉ riêng độ chính xác có thể đánh lừa.
• Để ý cảnh báo phân tách, đếm số sự kiện, và đừng khoác cho tỷ số odds vẻ ngoài của tỷ số nguy cơ.

Tiếp theo: cơ chế đằng sau mọi khoảng bạn đã thấy cho đến giờ, với khoảng tin cậy bằng t.test(), confint() và prop.test().