Exploring datasaurus

1 Giới thiệu

Như câu chuyện ngụ ngôn “Thầy bói sờ voi”, nhiều tác giả đã dùng những thí nghiệm mô phỏng để cảnh báo về nguy cơ sai lầm và nhược điểm của các phương pháp thống kê quy ước, bao gồm Mean, Sd, Pearson’s r, Boxplot..

Thí dụ:

Nghịch lý Simpson E.H (1951) là một can thiệp nhắm vào hệ số tương quan r của Pearson; theo đó : Trong dữ liệu được tái cấu trúc như hình B: Rõ ràng là Y tương quan nghịch với X, nhưng kết quả thống kê (Hình C) cho ra giá trị r>0 và tương đương với một dữ liệu có tương quan thuận trong hình A.

Năm 1973, ta biết đến bộ tứ Anscombe với 4 dữ liệu với cấu trúc hoàn toàn khác nhau nhưng đều cho ra cùng kết quả Trung bình, cùng độ lệch chuẩn và cùng hệ số tương quan. Dĩ nhiên ta có thể phủ nhận ý tưởng của Anscombe vì cho rằng nó quá cực đoan. Nhưng vào năm 2016, Justin Matekja và George Fitzmaurice (https://www.autodeskresearch.com/publications/samestats) đã sử dụng một thuật toán mô phỏng để hiện thực hóa bộ tứ Anscombe theo hướng ngẫu nhiên, phi cấu trúc. Kết quả cho thấy cả 8 bộ dữ liệu này đều có chung kết quả Mean, Sd, Pearson’s r. Kết quả này cho thấy là nguy cơ nhầm lẫn hoàn toàn có thể xảy ra trên thực tế. Một cách giải thích cho nghịch lý này, đó là “Regression towards the means” hay hiện tượng Hồi quy về trung bình.

2 12 biến thể Datasaurus

Cũng trong năm 2016, Cairo A. còn đi xa hơn nữa khi giới thiệu không chỉ 4, mà đến 12 biến thể dữ liệu. Thí nghiệm nổi tiếng này có tên là Datasaurus Dozen. Trong bài này, Nhi sẽ cùng các bạn khảo sát bộ dữ liệu này.

Trước hết ta tải file Datasaurus từ đây: https://www.autodeskresearch.com/sites/default/files/The%20Datasaurus%20Dozen.zip

library(tidyverse)

df=read_tsv("DatasaurusDozen.tsv")%>%as_tibble()

3 Các trị số thống kê mô tả

Bộ số liệu có dạng “long format”, gồm 3 biến, biến dataset cho phép phân chia dataframe thành 13 nhóm khác nhau, mỗi bộ dữ liệu nhỏ này có 2 biến liên tục X và Y, với số lượng như nhau (n=142).

Thống kê mô tả với Mean và SD cho thấy cả 13 bộ số liệu đều có Mean X = 54.27 và SD của X = 16.77

df %>%
  group_by(dataset) %>%
  summarise_at("x",funs(Size=length,
                        Mean=mean,
                        SD=sd))%>%knitr::kable()

dataset	Size	Mean	SD
away	142	54.26610	16.76983
bullseye	142	54.26873	16.76924
circle	142	54.26732	16.76001
dino	142	54.26327	16.76514
dots	142	54.26030	16.76774
h_lines	142	54.26144	16.76590
high_lines	142	54.26881	16.76670
slant_down	142	54.26785	16.76676
slant_up	142	54.26588	16.76885
star	142	54.26734	16.76896
v_lines	142	54.26993	16.76996
wide_lines	142	54.26692	16.77000
x_shape	142	54.26015	16.76996

Làm tương tự cho biến Y, ta sẽ thấy kết quả là Mean Y=47.83 và SD Y=26.94, như nhau cho 13 nhóm.

df %>%
  group_by(dataset) %>%
  summarise_at("y",funs(Size=length,
                        Mean=mean,
                        SD=sd))%>%knitr::kable()

dataset	Size	Mean	SD
away	142	47.83472	26.93974
bullseye	142	47.83082	26.93573
circle	142	47.83772	26.93004
dino	142	47.83225	26.93540
dots	142	47.83983	26.93019
h_lines	142	47.83025	26.93988
high_lines	142	47.83545	26.94000
slant_down	142	47.83590	26.93610
slant_up	142	47.83150	26.93861
star	142	47.83955	26.93027
v_lines	142	47.83699	26.93768
wide_lines	142	47.83160	26.93790
x_shape	142	47.83972	26.93000

Phân tích tương quan giữa Y và X cho thấy cả 13 nhóm đều có hệ số tương quan Pearson’s r rất gần nhau, khoảng từ -0.063 đến -0.069

df%>%split(.$dataset)%>%map(~cor(.$x,.$y,method="pearson"))%>%cbind()

##            .          
## away       -0.06412835
## bullseye   -0.06858639
## circle     -0.06834336
## dino       -0.06447185
## dots       -0.06034144
## h_lines    -0.06171484
## high_lines -0.06850422
## slant_down -0.06897974
## slant_up   -0.06860921
## star       -0.0629611 
## v_lines    -0.06944557
## wide_lines -0.06657523
## x_shape    -0.06558334

Xa hơn nữa, ta có thể thấy là 13 mô hình hồi quy Y~X đều có nội dung như nhau:

df%>%split(.$dataset)%>%
  map(~broom::tidy(lm(.$y~.$x)))

## $away
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.4251300 7.6932059  6.9444560 1.308597e-10
## 2         .$x -0.1030184 0.1354896 -0.7603419 4.483288e-01
## 
## $bullseye
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.8094708 7.6903593  6.9970035 9.917561e-11
## 2         .$x -0.1101675 0.1354339 -0.8134407 4.173467e-01
## 
## $circle
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.7970450 7.6925470  6.9933983 1.010830e-10
## 2         .$x -0.1098143 0.1354821 -0.8105447 4.190029e-01
## 
## $dino
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.4529784 7.6933924  6.9479075 1.285021e-10
## 2         .$x -0.1035825 0.1355026 -0.7644316 4.458966e-01
## 
## $dots
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.0983419 7.6924229  6.9026811 1.630167e-10
## 2         .$x -0.0969127 0.1354905 -0.7152729 4.756316e-01
## 
## $h_lines
##          term    estimate std.error  statistic      p.value
## 1 (Intercept) 53.21108724 7.6954597  6.9146080 1.531137e-10
## 2         .$x -0.09916499 0.1355427 -0.7316144 4.656268e-01
## 
## $high_lines
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.8087932 7.6926945  6.9947914 1.003417e-10
## 2         .$x -0.1100695 0.1354766 -0.8124615 4.179063e-01
## 
## $slant_down
##          term   estimate std.error statistic      p.value
## 1 (Intercept) 53.8497077 7.6911834  7.001485 9.685386e-11
## 2         .$x -0.1108172 0.1354522 -0.818128 4.146744e-01
## 
## $slant_up
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.8125956 7.6909632  6.9968604 9.925065e-11
## 2         .$x -0.1102184 0.1354513 -0.8137125 4.171915e-01
## 
## $star
##          term  estimate std.error  statistic      p.value
## 1 (Intercept) 53.326679 7.6915981  6.9331078 1.389166e-10
## 2         .$x -0.101113 0.1354591 -0.7464467 4.566492e-01
## 
## $v_lines
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.8908434 7.6903136  7.0076262 9.375965e-11
## 2         .$x -0.1115508 0.1354299 -0.8236796 4.115226e-01
## 
## $wide_lines
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.6349489 7.6914773  6.9732961 1.124023e-10
## 2         .$x -0.1069408 0.1354572 -0.7894803 4.311664e-01
## 
## $x_shape
##          term   estimate std.error  statistic      p.value
## 1 (Intercept) 53.5542263 7.6888696  6.9651625 1.173303e-10
## 2         .$x -0.1053169 0.1354267 -0.7776668 4.380777e-01

4 Sự thật về Datasaurus

Sự thật về cấu trúc của 13 bộ số liệu này chỉ được phơi bày khi ta mô tả trực quan với biểu đồ tán xạ. Bạn sẽ bất ngờ vì 13 hình ảnh hoàn toàn khác biệt:

df%>%ggplot(aes(x,y))+
  geom_point(aes(col=dataset),size=1)+
  theme_bw(8)+
  coord_equal()+
  facet_wrap(~dataset,ncol=4,scales="free")

Như vậy, chắc các bạn cũng đồng ý rằng những trị số Trung bình và SD hoàn toàn bị hack (đánh lừa), chúng không đủ để mô tả trung thực về đặc tính phân bố của dữ liệu.

Câu hỏi đặt ra đó là: Nếu Mean và SD có nguy cơ sai lầm cao như vậy, liệu ta có thể tránh nhưng ảo giác này bằng những phép thống kê mô tả khác hay không ?

Để trả lời, Nhi sẽ lần lượt thử nghiệm những cách làm khác nhau:

1. Đầu tiên, ta sử dụng thêm Median, IQR, Skewness, kurtosis cho biến Y

Kết quả cho thấy Median, Skewness, Kurtosis và IQR không bị ảnh hưởng bởi thủ thuật can thiệp trong 13 bộ dữ liệu datasaurus.

library(e1071)

df %>%
  group_by(dataset) %>%
  summarise_at("y",funs(Mean=mean,
                        Median=median,
                        Skewness=skewness,
                        Kurtosis=kurtosis,
                        IQR=IQR))%>%knitr::kable()

dataset	Mean	Median	Skewness	Kurtosis	IQR
away	47.83472	47.53527	-0.0932474	-1.2817225	47.17726
bullseye	47.83082	47.38294	0.0122276	-1.4540983	46.28812
circle	47.83772	51.02502	0.0048960	-1.8029697	59.43277
dino	47.83225	46.02560	0.2472603	-1.0635521	43.23723
dots	47.83983	51.29929	0.1476910	-1.4088664	65.77445
h_lines	47.83025	50.47353	0.1487355	-1.1603030	39.86956
high_lines	47.83545	32.49920	0.2024742	-1.8022348	53.01918
slant_down	47.83590	46.40131	0.1873369	-0.9492562	40.59857
slant_up	47.83150	45.29224	0.2554495	-1.0933159	46.09960
star	47.83955	50.11055	0.2372490	-1.3589160	43.17445
v_lines	47.83699	47.11362	0.2253746	-1.0479510	43.09251
wide_lines	47.83160	46.27933	0.2484350	-1.0901205	43.22119
x_shape	47.83972	39.87621	0.2277144	-1.4442001	50.13882

2. Bách phân vị, Min và Max

Tương tự, khảo sát bách phân vị cũng cho phép mô tả trung thực nhất đặc tính phân bố của biến Y

df %>%
  group_by(dataset) %>%
  summarise_at("x",funs(Min=min,Max=max,
                        Q5=quantile(.,probs=0.05),
                        Q25=quantile(.,probs=0.55),
                        Median=quantile(.,probs=0.5),
                        Q75=quantile(.,probs=0.75),
                        Q95=quantile(.,probs=0.95)))%>%knitr::kable()

dataset	Min	Max	Q5	Q25	Median	Q75	Q95
away	15.56075	91.63996	30.41458	59.84299	53.34030	69.14660	78.04134
bullseye	19.28820	91.73554	25.18011	55.11872	53.84209	64.79890	88.42245
circle	21.86358	85.66476	26.31463	56.32046	54.02321	64.97267	84.61663
dino	22.30770	98.20510	29.50002	55.26922	53.33330	64.74360	85.03845
dots	25.44353	77.95444	25.98733	51.88621	50.97677	75.19736	77.24589
h_lines	22.00371	98.28812	27.17405	54.89325	53.06968	66.76827	80.39223
high_lines	17.89350	96.08052	27.31751	56.54162	54.16869	63.95267	84.88861
slant_down	18.10947	95.59342	27.74931	55.07874	53.13516	64.46999	82.77544
slant_up	20.20978	95.26053	27.89781	55.57817	54.26135	64.48801	85.79441
star	27.02460	86.43590	28.59477	57.58533	56.53473	68.71149	79.16531
v_lines	30.44965	89.50485	30.47846	50.42004	50.36289	69.50407	89.49400
wide_lines	27.43963	77.91587	30.80572	65.29557	64.55023	67.45367	74.16283
x_shape	31.10687	85.44619	34.81317	50.75597	47.13646	71.85692	81.52456

3. Bootstrap trung bình

Ta cũng có thể áp dụng một phân tích trung bình kèm theo bootstrap. Việc chọn mẫu ngẫu nhiên có thể khắc phục nguy cơ sai lầm: Giá trị Mean tính được cho y sau 100 lần bootstrap không còn đồng nhất như ban đầu:

library(broom)

set.seed(1234)

bootmean=function(data,n){
  data%>%bootstrap(n)%>%
  do(tidy(lm(.$y~1)))->temp
  mean(temp$estimate)
}

df%>%split(.$dataset)%>%map(~bootmean(.,n=100))%>%cbind()

##            .       
## away       47.64329
## bullseye   47.96658
## circle     48.24137
## dino       47.47619
## dots       47.75413
## h_lines    47.98973
## high_lines 47.85328
## slant_down 48.03726
## slant_up   48.23992
## star       48.07929
## v_lines    47.96137
## wide_lines 47.75999
## x_shape    48.06804

4. Tương quan phi tham số:

Vì dữ liệu X và Y đều có phân bố bất thường, do đó khi ta sử dụng phương pháp tương quan Spearman’s rho, kết quả sẽ khác hẳn so với Pearson’s r.

df%>%split(.$dataset)%>%map(~cor(.$x,.$y,method="spearman"))%>%cbind()

##            .           
## away       -0.05729991 
## bullseye   -0.07873367 
## circle     -0.0772919  
## dino       -0.06510904 
## dots       -0.1263792  
## h_lines    -0.0519729  
## high_lines -0.002868872
## slant_down -0.06693546 
## slant_up   -0.0860976  
## star       -0.05144481 
## v_lines    -0.05662093 
## wide_lines -0.05223275 
## x_shape    -0.02053475

5 Mô tả trực quan

Kết luận hiển nhiên từ 2D scatter plot của Datasaurus cho thấy việc mô tả trực quan có vai trò rất quan trọng, số liệu thống kê là hư ảo, nhưng biểu đồ luôn nói lên sự thật. Tuy nhiên, Nhi muốn tìm hiểu xa hơn một chút, với câu hỏi: Ngoài 2D plot, khi muốn mô tả đặc tính phân bố 1 chiều (đơn biến), dạng biểu đồ nào là tốt nhất ?

Ta sẽ lần lượt thử các giải pháp sau đây cho biến Y:

1. Point range (Điểm và error bar)

Kết quả cho thấy Point range plot là dạng biểu dồ đơn giản nhất, nhưng nó hoàn toàn thất bại và bị đánh lừa :

df%>%ggplot(aes(x=dataset,y))+
  stat_summary(aes(col=dataset),
               fun.data = "mean_cl_boot",
               size=1)+theme_bw()+coord_flip()

2. Boxplot:

Boxplot là dạng biểu đồ rất thông dụng để mô tả đặc tính phân bố cho 1 biến. Tuy nhiên, chính nguyên lý đơn giản (chỉ sử dụng 5 thông tin) trong Boxplot khiến nó có nguy cơ bị ngộ nhận cao. Khi áp dụng cho 13 bộ dữ liệu Datasaurus, khả năng phân biệt của Boxplot rất kém, hơn nữa nó còn gây ảo giác về phân bố chuẩn, liên tục của biến Y trong một số trường hợp như nhóm h_lines, dots, circle, dino.

p1=df%>%ggplot(aes(x=dataset,y=y))+
  geom_boxplot(aes(fill=dataset),alpha=0.5,show.legend = F)+
  coord_flip()+
  theme_bw()

p2=df%>%ggplot(aes(x=dataset,y=y))+
  geom_jitter(aes(col=dataset),size=1,alpha=0.6,show.legend = F)+
  scale_x_discrete(breaks=NULL)+
  coord_flip()+
  theme_bw()

gridExtra::grid.arrange(p1,p2,ncol=2)

3. Violin plot:

Violin plot có ưu thế cao hơn hẳn so với boxplot, trong trường hợp này nếu sử dụng violin plot là có thể phân biệt dễ dàng đặc tính phân phối của Y giữa 13 bộ số liệu. Nhược điểm duy nhất của violin plot là nó vẫn còn bị đánh lừa bởi cấu trúc đối xứng trong các nhóm hlines và bullseye; kết quả sẽ gây ảo giác về một phân phối uniform trong khi sự thật không phải như vậy.

p1=df%>%ggplot(aes(x=dataset,y))+
  geom_violin(aes(fill=dataset),alpha=0.5,show.legend = F)+
  coord_flip()+
  theme_bw()

p2=df%>%ggplot(aes(x=dataset,y))+
  geom_jitter(aes(col=dataset),size=1,alpha=0.6,show.legend = F)+
  scale_x_discrete(breaks=NULL)+
  coord_flip()+
  theme_bw()

gridExtra::grid.arrange(p1,p2,ncol=2)

4. Histogram

Histogram là dạng biểu đồ tốt nhất để mô tả đặc tính phân phối đơn biến. Thật vậy, histogram cung cấp hình ảnh trung thực nhất, hơn nữa - có thể tinh chỉnh tùy ý - về mật độ phân bố.Ta có thể phân biệt dễ dàng và trực quan 13 phân nhóm.

library(ggridges)

p1=df%>%ggplot(aes(y=dataset,x=y))+
  geom_density_ridges(aes(fill=dataset),
                      alpha=0.5,
                      scale=4,
                      stat = "binline",
                      bins=50,
                      show.legend = F)+
  theme_bw()

p2=df%>%ggplot(aes(x=dataset,y))+
  geom_jitter(aes(col=dataset),size=1,alpha=0.6,show.legend = F)+
  scale_x_discrete(breaks=NULL)+
  coord_flip()+
  theme_bw()

gridExtra::grid.arrange(p1,p2,ncol=2)

5. Kernel density plot

Khi làm mượt histogram, ta sẽ có Density plot, tuy nhiên cũng chính vì lý do này, density plot dễ bị nhiễu hơn histogram. Cũng như violin plot, densityplot hoàn toàn bị đánh lừa với các phân phối đối xứng, multimodal và gây ảo giác về phân phối chuẩn. Do đó, densityplot kém hơn histogram.

library(ggridges)

p1=df%>%ggplot(aes(y=dataset,x=x))+
  geom_density_ridges(aes(fill=dataset),alpha=0.5,scale=3,show.legend = F)+
  theme_bw()

p2=df%>%ggplot(aes(x=dataset,y))+
  geom_jitter(aes(col=dataset),size=1,alpha=0.6,show.legend = F)+
  scale_x_discrete(breaks=NULL)+
  coord_flip()+
  theme_bw()

gridExtra::grid.arrange(p1,p2,ncol=2)

## Picking joint bandwidth of 5.46

6 Kết luận

Sau thí nghiệm nhỏ này, ta có thể rút ra một số kết luận như sau:

1. Bộ dữ liệu datasaurus chỉ là một can thiệp mô phỏng, tận dụng tối đa hiện tượng “Hồi quy về trung bình” và thủ thuật gây nhiễu để đánh lừa mô hình Hồi quy Tuyến tính và các công thức Mean, SD cho từng biến số. Dù cách làm này có vẻ cực đoan và phi thực tế, tuy nhiên nó đã chứng minh được rằng các bảng kết quả thống kê mô tả như Mean, SD mà ta thường thấy trong các công bố khoa học đều đáng ngờ.

2. Thống kê mô tả đơn giản bằng Mean và SD là không đủ để nói lên sự thật về cấu trúc dữ liệu. Để tránh nguy cơ bị đánh lừa bởi ảo giác, chúng ta cần phải khảo sát chi tiết hơn về đặc tính phân bố, thí dụ dùng thêm Min, Max, Skewness, Kurtosis, IQR và bách phân vị. Để an toàn hơn nữa, ta có thể sử dụng bootstrap.

3. Trong mọi trường hợp, thống kê mô tả luôn phải kết hợp với biểu đồ trực quan. Data visualisation rất quan trọng, thậm chí quan trọng hơn các bảng thống kê, vì chỉ có biểu đồ mới cho phép chúng ta nhìn thấy sự thật.

4. Cho mục tiêu mô tả đặc tính phân bố, ta có thể sử dụng nhiều giải pháp (dạng biểu đồ), tuy nhiên khả năng của mỗi dạng có thể khác nhau. Boxplot là dạng biểu đồ đơn giản và thông dụng nhất, nhưng nó cũng chính là phương pháp kém chính xác nhất. Violin plot chính xác và tinh tế hơn Boxplot. Histogram và scatter plot là 2 dạng biểu đồ chính xác nhất, gần với sự thực nhất. Density plot đẹp và giàu thông tin hơn boxplot, nhưng kém chính xác hơn Histogram.

Bài thực hành đến đây là kết thúc. Cảm ơn sự theo dõi của các bạn và hẹn gặp lại.

Tài liệu tham khảo:

Toàn bộ dataset có thể download tại đây:

https://www.autodeskresearch.com/sites/default/files/SameStatsDataAndImages.zip

1. Cairo, A. Download the Datasaurus: Never trust summary statistics alone; always visualize your data. http://www.thefunctionalart.com/2016/08/downloaddatasaurus-never-trust-summary.html.

2. Anscombe, F.J. (1973). Graphs in Statistical Analysis. The American Statistician 27, 1, 17–21.

3. Blyth, C.R. (1972). On Simpson’s Paradox and the Sure-Thing Principle. Journal of the American Statistical Association 67, 338, 364–366.

4. Justin Matejka, George Fitzmaurice. Same Stats, Different Graphs: Generating Datasets with Varied Appearance and Identical Statistics through Simulated Annealing. https://www.autodeskresearch.com/publications/samestats