Hồi quy tuyến tính trong R: lm(), summary() và predict()

Ý tưởng: đường bình phương nhỏ nhất

Hồi quy tuyến tính khớp một đường thẳng đi qua một đám điểm: y = intercept + slope × x. Trong mọi đường thẳng có thể, lm() chọn đường làm tổng bình phương phần dư nhỏ nhất; phần dư là khoảng cách theo chiều dọc giữa một điểm và đường thẳng. Bình phương khiến những lần lệch lớn bị tính nặng hơn hẳn, vì vậy một điểm ngoại lai cực đoan có thể làm nghiêng cả đường khớp.

Nếu tương quan cho bạn một con số không đơn vị cho câu hỏi "hai biến này đi cùng nhau chặt đến mức nào", thì hồi quy cho bạn một phương trình: có đơn vị, có độ dốc diễn giải được, và có công cụ để dự đoán.

fit <- lm(mpg ~ wt, data = mtcars)
print(coef(fit))

Đọc công thức là "mô hình hóa mpg như một hàm của wt". Hai hệ số chính là đường đã khớp: mpg ≈ 37,3 − 5,3 × cân nặng. Độ dốc có đơn vị thật: mỗi 1000 lbs cân nặng xe tăng thêm (wt tính theo đơn vị 1000 lbs) làm mất khoảng 5,3 dặm mỗi gallon. Hệ số chặn (37,3 mpg khi cân nặng bằng 0) chỉ là chỗ đường thẳng cắt trục tại 0; không có chiếc xe nào nặng 0, nên đừng diễn giải nó quá mức.

Đi qua từng phần của summary()

summary(fit) là đầu ra mà mọi khóa học thống kê đều yêu cầu bạn diễn giải. Chạy nó, rồi đọc từng khối một:

fit <- lm(mpg ~ wt, data = mtcars)
summary(fit)

Call: nhắc lại mô hình bạn đã khớp. Giờ thì tầm thường, nhưng sẽ cứu bạn khi phải xoay xở với sáu đối tượng mô hình cùng lúc.

Residuals: bản tóm tắt năm số của phần còn lại (thực tế − dự đoán). Bạn muốn trung vị gần 0 và Min/Max cũng như 1Q/3Q tương đối đối xứng; bất đối xứng mạnh gợi ý rằng mô hình đường thẳng đang bỏ sót điều gì đó.

Coefficients: trái tim của đầu ra, mỗi số hạng một hàng:

• Estimate: giá trị đã khớp. Với wt là −5,34: mỗi 1000 lbs tăng thêm gắn với việc giảm khoảng 5,3 mpg. Luôn chuyển độ dốc thành một câu có đơn vị; câu đó là toàn bộ nội dung thực tiễn của mô hình.
• Std. Error: ước lượng sẽ dao động bao nhiêu qua các mẫu lặp lại. Ước lượng chỉ cách 0 vài sai số chuẩn thì không vững.
• t value: Estimate ÷ Std. Error, tức hệ số nằm cách 0 bao nhiêu sai số chuẩn (ở đây là −9,56).
• Pr(>|t|): p-value cho câu hỏi "liệu hệ số này có thật sự bằng 0 không?". Với wt là khoảng 1,3e-10: nếu cân nặng thật sự không có quan hệ tuyến tính với mpg, một độ dốc dốc như vậy gần như không bao giờ xuất hiện trong một mẫu 32 xe. P-value nhỏ = bằng chứng rằng mối liên hệ tồn tại, không phải bằng chứng mô hình đúng, và cũng không phải thước đo mức độ quan trọng (một hiệu ứng rất nhỏ nhưng được ước lượng chính xác cũng có p-value rất nhỏ).
• Signif. codes / dấu sao: cách viết tắt trực quan cho cột p-value. Tiện lợi, nhưng không thêm thông tin gì.

Residual standard error: 3.05 on 30 degrees of freedom: độ lớn điển hình của một lần dự đoán lệch, tính bằng chính đơn vị của biến phản hồi: dự đoán thường lệch khoảng 3 mpg. Hãy đánh giá nó so với thang đo của mpg (dao động từ khoảng 10 đến 34).

Multiple R-squared: 0.75: cân nặng giải thích khoảng 75% phương sai của mpg. Adjusted R-squared (0.74) tính lại con số đó với một khoản phạt cho mỗi biến dự báo, vì phiên bản thô chỉ có thể tăng khi bạn thêm biến, kể cả nhiễu ngẫu nhiên. Khi so sánh các mô hình có số biến dự báo khác nhau, phiên bản hiệu chỉnh mới là con số trung thực. Và hãy cưỡng lại phản xạ "mô hình tốt = R² cao": một hiệu ứng thật sự hữu ích có thể nằm trong mô hình R² thấp (biến kết quả nhiều nhiễu, một trong nhiều yếu tố tác động), còn R² cao có thể đến từ quá khớp hoặc một biến bị rò rỉ thông tin.

F-statistic: 91.4 ... p-value: 1.29e-10: kiểm định toàn bộ mô hình: mô hình này có tốt hơn "cứ dự đoán giá trị trung bình cho tất cả" không? Với một biến dự báo, nó trùng với kiểm định t của độ dốc (để ý 9,56² ≈ 91,4); với nhiều biến dự báo, nó trở thành kiểm định đồng thời rằng ít nhất một hệ số khác 0. Cơ chế của nó là cùng một phép phân rã phương sai như ANOVA.

Hồi quy bội: giữ nguyên các biến khác

Thêm biến dự báo bằng +:

fit2 <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit2)

Cách diễn giải thay đổi ở một điểm then chốt. Mỗi Estimate giờ là tác động của biến dự báo đó khi giữ nguyên các biến còn lại: hệ số wt (khoảng −3,9, giảm từ −5,3) là mức mpg mất đi do thêm cân nặng khi so sánh các xe có cùng mã lực. Con số −5,3 của hồi quy đơn đã lặng lẽ gộp cả thực tế rằng xe nặng hơn cũng thường mạnh hơn; hồi quy bội tách điều đó ra. Đây cũng là lý do các hệ số thay đổi khi bạn thêm biến: nếu biến dự báo mới tương quan với một biến cũ, vai trò của biến cũ thay đổi theo. R bình phương tăng lên khoảng 0,83, và ở đây phiên bản hiệu chỉnh mới là phép so sánh công bằng với mô hình một biến dự báo.

Dự đoán: predict()

Mô hình đã khớp là một hàm; predict() tính giá trị của nó. Tạo một data frame newdata có tên cột khớp chính xác với các biến dự báo:

fit <- lm(mpg ~ wt, data = mtcars)

new_cars <- data.frame(wt = c(2.5, 3.5))
print(predict(fit, newdata = new_cars))

print(predict(fit, newdata = new_cars, interval = "confidence"))
print(predict(fit, newdata = new_cars, interval = "prediction"))

Hai loại khoảng trả lời hai câu hỏi khác nhau, và nhầm lẫn chúng là lỗi kinh điển trong bài thi:

• interval = "confidence": độ bất định về giá trị trung bình: "với mọi xe nặng 2500 lbs, mpg trung bình nằm ở đâu?". Hẹp, và thu hẹp dần khi dữ liệu nhiều lên.
• interval = "prediction": khoảng mà một chiếc xe mới riêng lẻ có cân nặng đó nhiều khả năng rơi vào. Rộng hơn nhiều, vì mỗi chiếc xe mang độ phân tán riêng quanh đường thẳng, thứ phân tán mà dù có bao nhiêu dữ liệu cũng không lấy trung bình mà triệt tiêu được.

Báo cáo khoảng tin cậy khi câu hỏi là về một quan sát mới sẽ phóng đại độ chính xác của bạn rất nhiều.

Chẩn đoán và bẫy ngoại suy

summary() cho bạn biết mô hình ước lượng gì; các biểu đồ phần dư cho bạn biết có nên tin nó không. Trong một phiên tương tác:

par(mfrow = c(2, 2))
plot(fit)   # four diagnostic plots

Cần tìm gì: Residuals vs Fitted nên là một đám mây không có hình dạng; một đường cong nghĩa là quan hệ không thẳng, một hình phễu (độ phân tán tăng theo giá trị dự đoán) nghĩa là phương sai không đổi bị vi phạm, và sai số chuẩn của bạn sẽ sai. Các điểm trên Q-Q plot nên bám sát đường thẳng; đuôi nặng nghĩa là có điểm ngoại lai đang làm méo mô hình. Scale-Location là phép kiểm tra hình phễu một lần nữa. Residuals vs Leverage đánh dấu các điểm có ảnh hưởng lớn, những quan sát mà chỉ riêng chúng đã kéo lệch các hệ số (trong mtcars, những xe đặc biệt như Chrysler Imperial hay xuất hiện ở đây). Một biểu đồ phân tán nhanh của dữ liệu thô trước khi khớp mô hình sẽ phát hiện sớm phần lớn những điều này.

Cuối cùng, cái bẫy mà không phép chẩn đoán nào bắt được: ngoại suy. Mô hình học từ những chiếc xe nặng khoảng từ 1500 đến 5400 lbs. Đưa cho predict() một wt bằng 8, nó sẽ vui vẻ trả về mpg âm: phép toán kéo dài đường thẳng mãi mãi, nhưng bằng chứng dừng lại ở rìa dữ liệu. Chỉ dự đoán trong (hoặc gần) khoảng dữ liệu bạn đã dùng để khớp.

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

• fit <- lm(y ~ x, data = df) khớp đường bình phương nhỏ nhất; coef(fit) là phương trình, summary(fit) là báo cáo đầy đủ.
• Đọc Estimate thành câu có đơn vị; Pr(>|t|) hỏi "có thể bằng 0 không?", không hỏi "có quan trọng không?".
• Residual standard error là mức lệch điển hình theo đơn vị thật; adjusted R-squared là con số công bằng để so sánh mô hình.
• Trong hồi quy bội, mỗi hệ số mang nghĩa "khi giữ nguyên các biến còn lại", và các hệ số thay đổi khi có biến dự báo tương quan tham gia.
• predict(fit, newdata, interval = ...): "confidence" cho giá trị trung bình, "prediction" cho một trường hợp mới, loại rộng hơn.
• Kiểm tra plot(fit) để tìm đường cong, hình phễu và điểm có ảnh hưởng lớn; đừng bao giờ tin dự đoán nằm ngoài khoảng dữ liệu.

Tiếp theo: khi biến kết quả là có/không thay vì một con số, với hồi quy logistic bằng glm().