Overplotting

1 Giới thiệu

Đây là bài đầu tiên trong series “ggplot2 thực dụng”. Như chúng ta đều biết, ggplot2 là một công cụ đồ họa thống kê rất tuyệt vời trong R. Khác với những software thương mại, ggplot2 không cung cấp những dạng biểu đồ thống kê quy ước, nhưng nó cho phép người dùng tự do diễn đạt ý tưởng của mình thông qua “ngữ pháp đồ thị”.

Tuy nhiên, chính vì tinh thần tự do này, rất hiếm giáo trình ggplot2 thực sự hướng dẫn cho người học cách vận dụng ngữ pháp nói trên để đạt những mục tiêu thực dụng. Đa số giáo trình chỉ mới dừng lại ở việc giải thích cú pháp các functions ggplot2 hoặc bám vào những dạng biểu đồ kinh điển. Ngay cả khi thông thạo ngữ pháp ggplot2, người dùng vẫn bị chi phối bởi thói quen và cách nhìn vấn đề. Thí dụ, khi bạn gặp vấn đề khó, bạn vẫn phải tự suy nghĩ tìm ra giải pháp cho riêng mình.

Một trong những vấn đề này, đó là: Trình bày thông tin khi dataset có quá nhiều cases. Đây là vấn đề “Overplotting”, tức sự chồng lắp của các thành phần hình họa (điểm) khi chúng có tọa độ trùng nhau hay rất gần nhau. Do làm việc trong ngành kỹ thuật y sinh, Nhi thường phải xử lý những dataset với hàng trăm ngàn, thậm chí hàng triệu lines. Những giải pháp mà Nhi sắp chia sẻ sau đây có thể giải quyết hiện tượng chồng lắp và tóm tắt được thông tin về khuynh hướng của dữ liệu.

Trong thí dụ minh họa sau đây, Nhi tạo ra 1 dataset với 2 biến X,Y có phân bố chuẩn, với 99,999 cases.

library(tidyverse)

df=data_frame(X=rnorm(n=99999,mean=20,sd=50),
              Y=rnorm(n=99999,mean=10,sd=40)+X)%>%
  as_tibble()

Đây là một con số khá lớn, nên khi vẽ scatter plot trong ggplot2 biểu đồ sẽ rất rối.

Hình 1: Như bạn thấy: mật độ điểm quá dày và chúng chồng lên nhau một cách hỗn loạn:

p=ggplot(df,aes(X,Y))+ theme_bw()

p+geom_jitter(col="red")

Giải pháp thứ nhất, đó là can thiệp trên chính yếu tố hình họa (điểm), cụ thể là kích thước, độ tương phản, độ trong suốt.

Để tăng độ tương phản và trong suốt, ta có thể dùng những điểm rỗng, tức không có màu nền mà chỉ có màu viền:

Hình 2: Loại bỏ màu nền có thể phân biệt khá nhiều điểm, tuy nhiên với cỡ mẫu quá lớn, cách làm này chưa giải quyết được vấn đề chồng lắp:

p+geom_jitter(col="red",shape=1)

Hình 3:

Cách làm thứ 2, đó là thu nhỏ kích thước mỗi điểm chỉ còn 1 pixel. Biểu đồ sẽ thay đổi theo hướng tăng độ phân giải: TUy nhiên, sự chồng lắp vẫn còn;

p+geom_jitter(col="red",shape=".")

Hình 4: Tiếp theo, ta có thể thử can thiệp cùng lúc lên alpha (độ trong suốt của yếu tố hình họa) và loại bỏ màu nền:

Hiệu ứng của việc giảm alpha là giảm nhiễu, và giữ lại khuynh hướng chính của dữ liệu.

p+geom_jitter(col="red",alpha=0.03,shape=1)

Giải pháp can thiệp lên hình họa dừng lại ở đây, bây giờ ta sẽ thử một cách khác, đó là can thiệp lên yếu tố Thống kê. Cụ thể, ta muốn gom các giá trị gần nhau vào những nhóm đại diện, và biểu diễn tọa độ của các nhóm đại diện này.

Hình 5: geom_bin2d

library(viridis)

## Loading required package: viridisLite

p+geom_bin2d(bins=33,aes(col=..count..),show.legend = F)+
  scale_fill_viridis(begin=0.1,end=0.9,option="A",alpha=0.1)+
  scale_color_viridis(option="B",alpha=0.3)

Hình 6: Hoặc sử dụng geom_hex

p+geom_hex(aes(col=..count..),show.legend = F)+
  scale_fill_viridis(option="A",alpha=0.1,begin=0.1,end=0.9)+
scale_color_viridis(option="B",alpha=0.3)

Hai cách làm trên sẽ gom nhiều điểm lại thành một nhóm, đại diện bằng một ô vuông hay hình lục giác, và tô màu theo mật độ phân bố của các điểm

Ta có thể thay đổi hoàn toàn yếu tố thống kê, thí dụ thể hiện mật độ phân bố 2 chiều thay vì tọa độ giá trị:

Hình 7: 2D Kernel density plot

p+stat_density_2d(aes(fill = ..level..),geom = "polygon",col="black")+
  scale_fill_distiller(palette="Spectral")

Hình 8: Tương tự như trên, Ta có thể sử dụng các nhóm đại diện, nhằm cung cấp cả 2 thông tin về tọa độ và mật độ phân bố.

p+stat_density_2d(aes(fill = ..level..),geom = "polygon",alpha=0.3,show.legend = F)+ 
stat_density_2d(geom = "point", 
                  aes(size = ..density..,fill=..density..),
                  shape=21,n = 18, 
                  contour = FALSE,col="black")+
  scale_fill_distiller(palette="Spectral")

Hình 9: Có bản chất cũng là biểu đồ mật độ phân phối 2 chiều, nhưng có dạng 1 heatmap:

p+stat_density_2d(geom = "raster", 
                  aes(fill = ..density..), 
                  contour = FALSE)+
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  scale_fill_viridis(option="A",alpha=0.2)

Hình 10: Một heatmap khác đơn giản hơn:

p+stat_density_2d(geom = "raster", 
                  aes(fill = ..density..), 
                  contour = FALSE)+
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  scale_fill_gradient(low="white",high="red")

Bây giờ ta nhìn vấn đề theo một hướng khác: Giả định nếu ta chỉ quan tâm đến giá trị của 1 biến Y, ta có thể phân nhóm X và sử dụng boxplot, violin plot hay point range như sau:

Hình 11: Phân bố của Y theo X, sau khi phân nhóm X thành 10 đoạn

p+geom_boxplot(aes(group = cut_interval(X, 10),
                   fill=cut_interval(X, 10)),show.legend = F)+
  scale_fill_viridis(option="A",alpha=0.8,discrete = T)

Hình 12: Violin plot

p+geom_violin(aes(group = cut_interval(X, 10),
                   fill=cut_interval(X,10)),show.legend = F)+
  scale_fill_viridis(option="A",alpha=0.8,discrete = T)

Hình 13: Point và error bar

ggplot(df,aes(x=cut_interval(X,10),y=Y))+
               stat_summary(aes(col=cut_interval(X,10)),
                        geom = "pointrange",
                        fun.y = median,
                        fun.ymax = max,
                        fun.ymin=min,
                        size=1)+
                       theme_bw()+
  scale_color_viridis(option="A",discrete = T)+
  scale_x_discrete(breaks=NULL)

Hình 14: Kết hợp scatter dot plot và tóm tắt trung vị của Y, theo 20 phân nhóm của X

ggplot(df,aes(x=cut_interval(X,20),y=Y))+
  geom_jitter(aes(col=Y),alpha=0.05,shape=1)+
  stat_summary(fun.y = "median", 
               geom = "point",
               shape=21,
               size=3,
               fill="white")+
  theme_bw()+
  scale_color_viridis(option="A")+
  scale_x_discrete(breaks=NULL)

Hình 15: Nếu ta đồng thời cắt cả X và Y thành 30 phân nhóm bằng nhau thì kết quả sẽ ra sao ?

ggplot(df,aes(x=cut_interval(X,30),
              y=cut_interval(Y,30)))+
  geom_jitter(aes(col=X),alpha=0.5,shape=".")+
  theme_bw()+
  scale_color_viridis(option="A")+
  scale_x_discrete(breaks=NULL)+
  scale_y_discrete(breaks=NULL)

Hình 16: Cuối cùng, ta có thể vẽ density plot của Y theo từng phân nhóm của X

library(ggridges)

ggplot(df,aes(y=cut_interval(X,15),x=Y))+
  geom_density_ridges(aes(fill=cut_interval(X,15),
                          col=cut_interval(X,15)),
                          scale=4,show.legend = F)+
  theme_bw()+
  scale_y_discrete(breaks=NULL)+
  coord_flip()+
  scale_fill_viridis(option="A",alpha=0.5,discrete = T)+
  scale_color_viridis(option="C",alpha=0.5,discrete = T,begin=1,end=0)

## Picking joint bandwidth of 9.34

Tạm biệt các bạn, hy vọng bài viết giúp các bạn mua vui được vài trống canh.

Hẹn gặp lại các bạn lần tới.