Tương quan trong R: cor(), cor.test() và ma trận tương quan

Tương quan giữa hai biến: cor()

Tương quan (correlation) đặt câu hỏi: khi một biến tăng, biến kia có xu hướng cũng tăng (dương), giảm (âm), hay chẳng liên quan gì (gần bằng 0)? Trong R chỉ cần một lần gọi:

print(cor(mtcars$wt, mtcars$mpg))

Kết quả là khoảng −0,87: xe càng nặng càng tốn nhiên liệu, và mối quan hệ này mạnh. Con số đó là hệ số tương quan Pearson, thường được ký hiệu là r.

Đọc giá trị r

Hệ số luôn nằm giữa −1 và +1. Dấu cho biết chiều; độ lớn cho biết mức độ mạnh:

| |r| | Cách đọc thông thường | | --- | --- | | 0,0 đến 0,2 | không đáng kể | | 0,2 đến 0,4 | yếu | | 0,4 đến 0,6 | trung bình | | 0,6 đến 0,8 | mạnh | | 0,8 đến 1,0 | rất mạnh |

Hãy xem các mức này là gợi ý để bàn luận, không phải luật: trong vật lý r bằng 0,6 là đáng thất vọng, trong tâm lý học thì đó là thành tựu đáng tự hào. Hai tính chất nên ghi nhớ: r không có đơn vị (tương quan giữa cân nặng tính bằng tấn và mpg cho cùng r với cân nặng tính bằng kilogam và mpg, vì r được tính trên giá trị đã chuẩn hóa), và r chỉ đo mối liên hệ tuyến tính: một quan hệ hình chữ U hoàn hảo vẫn có thể có r ≈ 0.

Và câu nhất định phải nói: tương quan không phải là nhân quả. Doanh số kem và số ca đuối nước tương quan mạnh theo các tháng trong năm, không phải vì kem làm người ta chết đuối, mà vì mùa hè đẩy cả hai lên. Tương quan cho biết hai biến cùng biến động; nó không nói gì về lý do. Có thể x tác động lên y, có thể y tác động lên x, có thể một yếu tố thứ ba (ở đây là mùa) tác động lên cả hai. Để phân định cần đến thí nghiệm hoặc lập luận nhân quả cẩn thận, không phải một giá trị r lớn hơn.

Spearman và Kendall: khi Pearson không phải công cụ phù hợp

Pearson làm việc trên giá trị gốc, nên nó nhạy với giá trị ngoại lai và không thấy được các quan hệ cong. Đối số method chuyển sang các phương án dựa trên thứ hạng:

x <- c(1, 2, 3, 4, 5, 6, 7, 8)
y <- c(1, 4, 9, 16, 25, 36, 49, 300)  # increasing, with a wild last value

print(cor(x, y))                        # Pearson - dragged by the outlier
print(cor(x, y, method = "spearman"))   # ranks: perfectly monotonic
print(cor(x, y, method = "kendall"))

Spearman thay mỗi giá trị bằng thứ hạng của nó rồi tính Pearson trên các thứ hạng đó. Vì y luôn tăng khi x tăng, mọi thứ hạng đều khớp nhau và Spearman cho đúng 1: độ lớn của giá trị ngoại lai không còn quan trọng, chỉ vị trí của nó mới quan trọng. Hãy dùng Spearman khi dữ liệu có thứ bậc (thang đo khảo sát), lệch nhiều, hoặc khi quan hệ đơn điệu nhưng không thẳng. Kendall trả lời câu hỏi tương tự dựa trên các cặp cùng chiều và ngược chiều; nó ổn định hơn với mẫu nhỏ nhưng chậm hơn, và Spearman là lựa chọn mặc định phổ biến.

Ma trận tương quan

Để xem quan hệ giữa nhiều biến cùng lúc, hãy đưa cho cor() vài cột số:

m <- cor(mtcars[, c("mpg", "wt", "hp", "disp")])
print(round(m, 2))

Mỗi biến được so với mọi biến khác, với số 1 trên đường chéo (mọi biến đều tương quan hoàn hảo với chính nó) và hai nửa đối xứng như soi gương qua đường chéo. Làm tròn đến hai chữ số thập phân quan trọng hơn bạn nghĩ: ma trận chưa làm tròn là một bức tường chữ số, trong khi mục đích của ma trận là để lướt qua. Ở đây, lướt qua sẽ thấy mpg tương quan âm với cả ba biến kia (xe nặng hơn, mạnh hơn, động cơ lớn hơn thì tốn nhiên liệu hơn), còn wt, hp và disp đều tương quan dương mạnh với nhau: một cụm biến "xe to" sẽ trở nên quan trọng khi bạn tới hồi quy tuyến tính và cơn đau đầu đa cộng tuyến của nó.

Giá trị thiếu: đối số use

Khi có dữ liệu thiếu, mặc định của cor() là trả về NA thay vì đoán:

d <- mtcars[, c("mpg", "wt", "hp")]
d$wt[c(3, 7)] <- NA  # knock two values out

print(cor(d)[1, 2])
print(cor(d, use = "complete.obs")[1, 2])
print(cor(d, use = "pairwise.complete.obs")[1, 2])

• use = "complete.obs" bỏ mọi hàng chứa bất kỳ NA nào, rồi tính cả ma trận từ các hàng còn lại. Nhất quán nhưng lãng phí: một giá trị wt bị thiếu cũng loại hàng đó khỏi cặp mpg và hp.
• use = "pairwise.complete.obs" tính mỗi ô từ mọi hàng mà cặp đó có đủ dữ liệu. Cách này giữ được nhiều dữ liệu hơn, nhưng các ô khác nhau dựa trên các tập con khác nhau, nên đôi khi có thể tạo ra một ma trận không nhất quán nội tại.

Với vài giá trị NA lẻ tẻ thì cách nào cũng được; chỉ cần nói rõ bạn đã dùng cách nào.

Có ý nghĩa thống kê không? cor.test()

cor() cho một con số nhưng không cho biết liệu nó có thể chỉ là nhiễu hay không. cor.test() thực hiện kiểm định giả thuyết:

print(cor.test(mtcars$wt, mtcars$mpg))

Đi qua đầu ra từng khối một:

• t = −9.56, df = 30: thống kê kiểm định. Giả thuyết không là tương quan thật bằng 0; r quan sát được được chuyển thành thống kê t với n − 2 bậc tự do (32 xe − 2).
• p-value = 1.29e-10: nếu tương quan thật bằng 0, xác suất thấy một giá trị r xa 0 đến mức này trong mẫu 32 xe là khoảng 0,0000000001. Đó là bằng chứng áp đảo rằng mối liên hệ là có thật, dù hãy nhớ rằng p-value chỉ nói về việc r có khác 0 hay không, không nói mối quan hệ có lớn hay có nhân quả hay không.
• 95 percent confidence interval: −0.93 to −0.74: dải hợp lý cho tương quan thật. Thường hữu ích hơn p-value: ngay cả đầu lạc quan nhất của khoảng này vẫn là tương quan âm mạnh.
• sample estimates: cor = −0.87: chính con số mà cor() đã cho bạn.

Với mẫu nhỏ cần thận trọng hơn: khi n = 10, các tương quan ±0,5 xuất hiện do may rủi thường đến đáng lo, và khoảng tin cậy rộng sẽ cho bạn thấy điều đó. Hãy báo cáo khoảng tin cậy, không chỉ p-value.

Nhìn tận mắt: luôn vẽ biểu đồ

Hệ số tương quan nén cả một mối quan hệ thành một con số, và sự nén đó có thể che giấu độ cong, các cụm, hoặc một điểm duy nhất gánh hết mọi thứ. Trước khi tin bất kỳ giá trị r nào, hãy xem biểu đồ phân tán:

plot(mtcars$wt, mtcars$mpg)              # one pair
pairs(mtcars[, c("mpg", "wt", "hp", "disp")])  # every pair in the matrix

pairs() vẽ một lưới biểu đồ phân tán tương ứng với ma trận tương quan, cách nhanh nhất để kiểm tra rằng các con số có đúng ý nghĩa bạn nghĩ. Để có đồ họa ma trận kiểu heatmap chỉn chu, package corrplot (install.packages("corrplot"), rồi corrplot(cor(m))) là công cụ tiêu chuẩn.

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

• cor(x, y) cho r của Pearson: dấu là chiều, độ lớn là mức độ mạnh, luôn trong [−1, 1], không có đơn vị.
• Tương quan đo sự cùng biến động tuyến tính và không nói gì về nhân quả: luôn có thể có một biến thứ ba ẩn phía sau.
• method = "spearman" cho thứ hạng: dữ liệu thứ bậc, giá trị ngoại lai, quan hệ đơn điệu nhưng cong.
• cor(df) trên các cột số cho ma trận; hãy round(, 2) nó, và chú ý đối số use = khi có dữ liệu thiếu.
• cor.test(x, y) thêm p-value và khoảng tin cậy; hãy báo cáo khoảng tin cậy.
• Luôn xem biểu đồ phân tán trước khi tin vào con số.

Tiếp theo: khi câu hỏi chuyển từ "chúng có cùng biến động không?" sang "trung bình của nhóm này có khác nhóm kia không?", đó là kiểm định t.