ANSCOMBE’S QUARTET

Linear Regression Model: Bẫy Model và Tóm tắt Dữ Liệu

1 Bẫy Ý Nghĩa Thống Kê

Khi phân tích dữ liệu để xây dựng những linear regression models, nếu chúng ta tìm được một mô hình có ý nghĩa thống kê và đường biểu diễn y = a + bx phù hợp với phân bố dữ liệu, lúc đó chúng ta hài lòng với kết quả của mình và ít khi có bất kì hoài nghi nào về tính duy nhất đúng của chúng nữa.

Ví dụ kinh điển, ANSCOMBE’s QUARTET sẽ cho chúng ta biết là mình không nên chủ quan với những con số tóm tắt và ngưỡng ý nghĩa thống kê vừa đạt được mà hãy đặt vấn đề tiếp là liệu chúng ta đang bị vướng vào cái bẫy dữ liệu nào không?

Ví dụ được trình bày dưới đây sẽ giải thích cho chúng ta một vài điều thú vị.

2 Chuẩn bị dữ liệu

Tôi sử dụng dataset từ package MASS của ví dụ kinh điển ANSCOMBE’s QUARTET bao gồm 11 quan sát và 8 biến số x1:x4 và y1:y4.

3 Mô hình hồi qui tuyến tính y = a + bx

Trước hết chúng ta xem xét mối quan hệ giữa cặp biến số đầu tiên x1 và y1

library(ggplot2)
m1<-lm(y1~x1,data = anscombe)
summary(m1)

## 
## Call:
## lm(formula = y1 ~ x1, data = anscombe)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92127 -0.45577 -0.04136  0.70941  1.83882 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0001     1.1247   2.667  0.02573 * 
## x1            0.5001     0.1179   4.241  0.00217 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6665, Adjusted R-squared:  0.6295 
## F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217

p1 <- ggplot(anscombe, aes(x1, y1)) + 
  geom_point()+
 geom_smooth(method = lm, se = TRUE)
p1

Các thông số tương ứng với a = 3.0001 và b = 0.5001 cũng như mô hình hồi quy tuyến tính trên đều đạt ngưỡng ý nghĩa thống kê, Pr = 0.00217 với Adj. R-squared = 0.6295

Mô hình phân bố scarted giữa x1 và y1 nói lên mối quan hệ linear giữa hai biến số này vì chúng phân bố khá đồng nhất trong toàn dãi x1 và y1.

Đến đây, chúng ta thường khá hài lòng với kết quả của mình về mối quan hệ tuyến tính giữa hai biến số mà không hề đặt nghi vấn về những tình huống khác có thể xảy ra.

4 Những tình huống khác thường

Ngoài cặp biến số (x1,y1) ra, chúng ta tiếp tục phân tích mối quan hệ giữa các cặp biến số còn lại là (x2,y2), (x3,y3) và (x4,y4).

Trước hết xem means của chúng có khác nhau hay không.

summary(anscombe)

##        x1             x2             x3             x4    
##  Min.   : 4.0   Min.   : 4.0   Min.   : 4.0   Min.   : 8  
##  1st Qu.: 6.5   1st Qu.: 6.5   1st Qu.: 6.5   1st Qu.: 8  
##  Median : 9.0   Median : 9.0   Median : 9.0   Median : 8  
##  Mean   : 9.0   Mean   : 9.0   Mean   : 9.0   Mean   : 9  
##  3rd Qu.:11.5   3rd Qu.:11.5   3rd Qu.:11.5   3rd Qu.: 8  
##  Max.   :14.0   Max.   :14.0   Max.   :14.0   Max.   :19  
##        y1               y2              y3              y4        
##  Min.   : 4.260   Min.   :3.100   Min.   : 5.39   Min.   : 5.250  
##  1st Qu.: 6.315   1st Qu.:6.695   1st Qu.: 6.25   1st Qu.: 6.170  
##  Median : 7.580   Median :8.140   Median : 7.11   Median : 7.040  
##  Mean   : 7.501   Mean   :7.501   Mean   : 7.50   Mean   : 7.501  
##  3rd Qu.: 8.570   3rd Qu.:8.950   3rd Qu.: 7.98   3rd Qu.: 8.190  
##  Max.   :10.840   Max.   :9.260   Max.   :12.74   Max.   :12.500

Kết quả cho thấy các cặp (x,y) có n và mean khá giống nhau. Hay nói cách khác đây có thể xem là là 4 subsets được rút ra ngẫu nhiên từ cùng một dataset.

Chúng ta xem mối quan hệ của chúng thế nào giữa 3 cặp (x,y) còn lại.

m1<-lm(y1~x1,data=anscombe)
m2<-lm(y2~x2,data=anscombe)
m3<-lm(y3~x3,data=anscombe)
m4<-lm(y4~x4,data=anscombe)
summary(m1)

## 
## Call:
## lm(formula = y1 ~ x1, data = anscombe)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.92127 -0.45577 -0.04136  0.70941  1.83882 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0001     1.1247   2.667  0.02573 * 
## x1            0.5001     0.1179   4.241  0.00217 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6665, Adjusted R-squared:  0.6295 
## F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217

summary(m2)

## 
## Call:
## lm(formula = y2 ~ x2, data = anscombe)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9009 -0.7609  0.1291  0.9491  1.2691 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    3.001      1.125   2.667  0.02576 * 
## x2             0.500      0.118   4.239  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.237 on 9 degrees of freedom
## Multiple R-squared:  0.6662, Adjusted R-squared:  0.6292 
## F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179

summary(m3)

## 
## Call:
## lm(formula = y3 ~ x3, data = anscombe)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1586 -0.6146 -0.2303  0.1540  3.2411 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0025     1.1245   2.670  0.02562 * 
## x3            0.4997     0.1179   4.239  0.00218 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared:  0.6663, Adjusted R-squared:  0.6292 
## F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176

summary(m4)

## 
## Call:
## lm(formula = y4 ~ x4, data = anscombe)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.751 -0.831  0.000  0.809  1.839 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.0017     1.1239   2.671  0.02559 * 
## x4            0.4999     0.1178   4.243  0.00216 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.236 on 9 degrees of freedom
## Multiple R-squared:  0.6667, Adjusted R-squared:  0.6297 
## F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165

Nhìn vào summary của bốn models m1, m2, m3 và m4 ở các coefficients, Pr và Adj. R-squared ta thấy chúng hầu như tương đồng. Điều này cho chúng ta ý nghĩ là các cặp biên số (x,y) này đều có quan hệ linear với phương trình gần giống:

y = 3 + 0.5*x

Đây chính là cái bẫy dữ liệu.

Các bạn hãy nhìn vào biểu đồ phân bố và đường hồi qui của tất cả 4 cặp (x,y) khi được sắp đặt trên một panel nhé.

library(ggpubr)
library(gridExtra)
library(ggplot2)

p1 <- ggplot(anscombe) + geom_point(aes(x1, y1), color = "darkorange", size = 3) + theme_bw() + scale_x_continuous(breaks = seq(0, 20, 2)) + scale_y_continuous(breaks = seq(0, 12, 2)) + geom_abline(intercept = 3, slope = 0.5, color = "cornflowerblue") + expand_limits(x = 0, y = 0) + labs(title = "dataset 1")
p2 <- ggplot(anscombe) + geom_point(aes(x2, y2), color = "darkorange", size = 3) + theme_bw() + scale_x_continuous(breaks = seq(0, 20, 2)) + scale_y_continuous(breaks = seq(0, 12, 2)) + geom_abline(intercept = 3, slope = 0.5, color = "cornflowerblue") + expand_limits(x = 0, y = 0) + labs(title = "dataset 2")
p3 <- ggplot(anscombe) + geom_point(aes(x3, y3), color = "darkorange", size = 3) + theme_bw() + scale_x_continuous(breaks = seq(0, 20, 2)) + scale_y_continuous(breaks = seq(0, 12, 2)) + geom_abline(intercept = 3, slope = 0.5, color = "cornflowerblue") + expand_limits(x = 0, y = 0) + labs(title = "dataset 3")
p4 <- ggplot(anscombe) + geom_point(aes(x4, y4), color = "darkorange", size = 3) + theme_bw() + scale_x_continuous(breaks = seq(0, 20, 2)) + scale_y_continuous(breaks = seq(0, 12, 2)) + geom_abline(intercept = 3, slope = 0.5, color = "cornflowerblue") + expand_limits(x = 0, y = 0) + labs(title = "dataset 4")

grid.arrange(p1, p2, p3, p4, top = "Anscombe's Quartet")

Đến đây thì chúng ta đã hiểu vì sao nhà thống kê học Francis Anscombe lại tạo ra bộ dữ liệu này để cảnh báo chúng ta về cái bẫy ý nghĩa thống kê của các mô hình linear regression.

Trong cặp (x2, y2) và (x4, y4) thì mối quan hệ không thể gọi là linear nhưng mọi chỉ số đều cho thấy chúng quan hệ tuyến tính có ý nghĩa thống kê Pr và Adj. R-squared bảo đảm cả.

Thông điệp mà ông Anscombe muốn gởi đến chúng ta là hãy nhìn vào biểu đồ phân bố của (x,y) trước hết để có ý niệm ban đầu về mối quan hệ của chúng. Sau khi dựng được mô hình tuyến tính thì hãy đánh giá nó trong các datasets khác và so sánh các tiêu chí khác như variance, AIC… nhằm tránh được các loại bẫy ý nghĩa thống kê mà chúng ta mắc phải khi vô tình vi phạm một giả định nào đó của phép kiểm định.

5 Phân tích dữ liệu gộp lại

Để có cái nhìn tổng thể về toàn bộ các cặp (x,y) chúng ta gộp chung lại và phân tích như sau:

x<-c(anscombe$x1,anscombe$x2,anscombe$x3,anscombe$x3)
y<-c(anscombe$y1,anscombe$y2,anscombe$y3,anscombe$y4)
df<-data.frame(x,y)
dim(df)

## [1] 44  2

lm(y~x,data=df)

## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Coefficients:
## (Intercept)            x  
##      4.5587       0.3269

library(ggplot2)
p <- ggplot(df, aes(x, y)) + 
  geom_point()+
 geom_smooth(method = lm, se = TRUE)
p

Bạn nhận xét thế nào về phân bố của (x,y) khi gộp chung lại ? Khi có một dataset lớn hơn thì các mối quan hệ sẽ được chứng minh thuyết phục hơn.

Chân thành cám ơn các bạn đã đọc bài viết này. Bất cứ góp ý hay bàn luận nào của các bạn cũng luôn được mong đợi.