Thư viện cần thiết

Nếu chưa cài package pacman, bạn dùng lệnh install.packages("pacman")

pacman::p_load(tidyverse, ggplot2)

Trong đó:

• tidyverse: giúp biến đổ số liệu tốt hơn.
• ggplot2: giúp vẽ biểu đồ.

Cấu trúc cơ bản của ggplot2

Một câu lệnh đầy đủ để vẽ đượ một biểu đồ cơ bản trong ggplot2 là: ggplot(data, mapping) + geom_biểuđồ()

Trong đó:

• data: data frame chứa số liệu cần vẽ biểu đồ.
• mapping: đối tượng chứa các thông số về trục x, y, nhóm, màu sắc, kích thước,… được tạo ra từ aes().
• ggplot(data, mapping): nền cơ sở, từ đây sẽ dựng các biểu đồ khác lên.
• geom_biểuđồ(): loại biểu đồ cần vẽ. Ví dụ muốn vẽ biểu đồ tán xạ, ta dùng geom_point().

Có thể vẽ nhiều biểu đồ khác nhau cùng lúc trên cùng một nền cơ sở: ggplot(data, mapping) + geom_point + geom_bar +…

Ví dụ

# Số liệu sử dụng là mtcars có sẵn trong ggplot2
mtcars

##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
## Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

# Vẽ thử biểu đồ đầu tiên
ggplot(data = mtcars, mapping = aes(x = wt, y = mpg)) +
  geom_point(aes(size = qsec))

Nếu số liệu không được khai báo trong mỗi hàm vẽ biểu đồ (geom_point(), geom_bar(),…) thì hàm sẽ lấy số liệu từ nền cơ sở ggplot().

# Ví dụ vẽ nhiều biểu đồ trên cúng 1 nền cơ sở
ggplot(data = mtcars, mapping = aes(x = wt, y = mpg)) +
  geom_point() +
  geom_smooth(method="lm")

## `geom_smooth()` using formula 'y ~ x'

Lúc này geom_point() và geom_smooth() sẽ lấy số liệu của x, y trong ggplot(). Ta thấy có cả 2 loại biểu đồ trên cùng 1 nền cơ sở (cùng hệ trục tọa độ)

Tóm tắt toàn bộ lệnh trong ggplot2

Đây là bản tóm tắt (cheat sheet) mà mọi người có thể tải về, xem lúc cần thiết. Tải về tại đây.

Vẽ biểu đồ phân tán (scatter plot)

Biểu đồ phân tán (scatter plot):

• Sử dụng các dấu chấm để thể hiện giá trị (điểm giao nhau) của hai biến số khác nhau.
• Dùng để quan sát mối tương quan giữa 2 biến định lượng.
• Trục hoành (trục X) mô tả biến độc lập. Trục tung (Y) mô tả biến phụ thuộc.

ggplot(data=iris, mapping = aes(x = Sepal.Length, y = Sepal.Width)) +
  geom_point()

Lưu ý:

Không phải tất cả mối liên quan thể hiện trên biểu đồ phân tán đều là mối quan hệ nhân quả. Có thể do ảnh hưởng của các biến số khác khiến ta lầm tưởng sự tương quan khi khảo sát 2 biến số đơn độc.

head(mpg)

## # A tibble: 6 × 11
##   manufacturer model displ  year   cyl trans      drv     cty   hwy fl    class 
##   <chr>        <chr> <dbl> <int> <int> <chr>      <chr> <int> <int> <chr> <chr> 
## 1 audi         a4      1.8  1999     4 auto(l5)   f        18    29 p     compa…
## 2 audi         a4      1.8  1999     4 manual(m5) f        21    29 p     compa…
## 3 audi         a4      2    2008     4 manual(m6) f        20    31 p     compa…
## 4 audi         a4      2    2008     4 auto(av)   f        21    30 p     compa…
## 5 audi         a4      2.8  1999     6 auto(l5)   f        16    26 p     compa…
## 6 audi         a4      2.8  1999     6 manual(m5) f        18    26 p     compa…

mpg %>% 
  dplyr::group_by(class) %>%
  summarise(mean_cty = mean(cty), mean_hwy = mean(hwy)) %>% 
  ggplot(aes(x = mean_cty, y = mean_hwy)) +
  geom_point()

mpg %>% 
  dplyr::group_by(class) %>%
  summarise(mean_cty = mean(cty), mean_hwy = mean(hwy)) %>% 
  ggplot(aes(x = mean_cty, y = mean_hwy)) +
  geom_label(aes(label=class))

library(ggrepel)
mpg %>% 
  dplyr::group_by(class) %>%
  summarise(mean_cty = mean(cty), mean_hwy = mean(hwy)) %>% 
  ggplot(aes(x = mean_cty, y = mean_hwy)) +
  geom_label_repel(aes(label=class))

Vẽ biểu đồ histogram

Biểu đồ histogram dùng để mô tả phân bố của một biến liên tục (continuous).

ggplot(data = diamonds) +
  geom_histogram(aes(x = carat))

Ta chỉ cần cung cấp giá trị cho trục x, trục y sẽ được tự động tạo ra bằng các đếm số quan sát nằm trên từng khoảng của x (khoảng này là độ rộng của mỗi cột, còn gọi là binwidth)

# Thử binwidth = 1
ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = carat), binwidth = 1)

# Thử binwidth = 0.5
ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = carat), binwidth = 0.5)

ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = carat), bins = 4)

ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = price, fill = cut), binwidth = 1000)

Vẽ biểu đồ cột (bar chart)

Biểu đồ cột mô tả các nhóm trong một biến (categorical) theo một đại lượng khác, có thể là tần số của nhóm đó hoặc môt giá trị đo lường khác.

# Ví dụ 1
ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut))

Ví dụ 1 thể hiện số lượng của từng nhóm của biến cut. Ta dễ dàng thấy số lượng của các nhóm tăng dẩn từ trái qua phải. Nhóm ideal có số lượng nhiều nhất.

# Ví dụ 2
ggplot(data = pressure) +
  geom_col(mapping = aes(x = temperature, y = pressure))

Ví dụ 2 cho thấy áp suất đo được trên từng khoảng nhiệt độ khác nhau. Ta cũng dễ dàng nhân ra được xu hướng “Nhiệt độ (temperature) càng tăng thì áp suất (pressure) cũng tăng theo”.

2 ví dụ trên sử dụng 2 hàm khác nhau để vẽ biểu đồ cột là geom_bar() và geom_col(). geom_bar() không có tham số y vì trục y chính là số lượng của mỗi nhóm trên trục x.

# color
ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut, color = cut))

# fill
ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut, fill = cut))

# width
# Lưu ý: width không nằm trong aes() mà nằm ở ngoài nhé
ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut, fill = cut), width = 0.5)

ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut, fill = clarity))

ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut, fill = clarity), position = "dodge")

ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut, fill = clarity), position = "fill")

ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut, fill = clarity), position = "dodge") +
  coord_flip()

Vẽ biểu đồ quạt (pie chart)

Đây là biểu đồ không có ý nghĩa trong các báo cáo khoa học vì giá trị của nó không hơn gì một bảng số liệu nhỏ gọn, đơn giản. Ví vậy ggplot không hỗ trợ vẽ biểu đồ quạt. Tuy nhiên ta vẫn có mẹo để vẽ được. Bằng các thêm hàm coord_polar() vào trong biểu đồ cột.

diamonds %>% 
  select(cut) %>% 
  group_by(cut) %>% 
  count() %>% 
  ggplot(aes(x="", y=n, fill=cut)) +
  geom_bar(stat="identity", width=1) +
  coord_polar("y", start=0) + 
  theme_void() # Bỏ bớt các giá trị thừa

Vẽ biểu đồ hộp (boxplot)

Boxplot là biểu đồ dùng để mô tả một biến định lượng theo một cách trực quan và đơn giản. Qua boxplot, chúng ta sẽ có thông tin về:

• Median: Trung vị
• Quartile: Tứ phân vị
• IQR: Khoảng tứ phân vị
• Min, max: giá trị lớn nhất, giá trị nhỏ nhất (không kể các outliers)
• Outliers: Biến ngoại lai, giá trị nằm ngoài khoảng (Q1-1.5IQR; Q3+1.5IQR)

Ngoài ra ta có thể kiểm tra phân bố của biến liên tục có phải phân bố chuẩn hay không thông qua tính đối xứng trục của biểu đồ boxplot. Hình ảnh càng đối xứng thì phân bố của biến càng giống phân bố chuẩn.

ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = class, y = hwy))

Khi có nhiều nhóm, chúng ta có thể sử dụng boxplot nằm ngang để tiện hơn cho việc trình bày báo cáo. Bạn có thể làm điểu này bằng cách thay đổi x và y cho nhau, hoặc thêm hàm coord_flip() để xoay ngang biểu đồ.

# Cách 1
ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = hwy, y = class))

# Cách 2
ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = class, y = hwy)) +
  coord_flip()

Các biểu đồ tương tự boxplot: geom_dot(), geom_violin(). Bạn có thể xem thêm tại bài viết tham khảo.

Vẽ biểu đồ vẽ biểu đồ forrest

Biểu đồ forest thường được dùng trong các nghiên cứu tổng hợp meta-analysis. Nó giúp thể hiện tóm tắt kết quả của nhiều nghiên cứu khác nhau có cùng chung vấn đề nghiên cứu. Nếu các bạn chưa hiểu về forest plot có thể tham khảo 2 bài viết sau của R stats:

• Bài 1: Cách đọc biểu đồ forest
• Bài 2: Ứng dụng của biểu đồ Forest trong các nghiên cứu y học

# install.packages("devtools")
devtools::install_github("NightingaleHealth/ggforestplot")

df <-
  ggforestplot::df_linear_associations %>%
  filter(
    trait == "BMI",
    dplyr::row_number() <= 30
  )
df

## # A tibble: 30 × 5
##    name          trait    beta      se    pvalue
##    <chr>         <chr>   <dbl>   <dbl>     <dbl>
##  1 Isoleucine    BMI    0.339  0.00945 1.11e-281
##  2 Leucine       BMI    0.343  0.00951 1.25e-285
##  3 Valine        BMI    0.287  0.00951 7.94e-200
##  4 Phenylalanine BMI    0.343  0.00862 0        
##  5 Tyrosine      BMI    0.261  0.00900 6.65e-185
##  6 Alanine       BMI    0.179  0.00890 8.62e- 90
##  7 Glutamine     BMI   -0.134  0.00945 7.68e- 46
##  8 Glycine       BMI   -0.0296 0.00937 1.56e-  3
##  9 Histidine     BMI    0.0364 0.00917 7.25e-  5
## 10 Lactate       BMI    0.131  0.00911 9.20e- 47
## # … with 20 more rows

ggforestplot::forestplot(
  df = df,
  name = name,
  estimate = beta,
  se = se
)

# Thêm mức ý nghĩa, tên trục x và tên biểu đồ
ggforestplot::forestplot(
  df = df,
  estimate = beta,
  pvalue = pvalue,
  psignif = 0.002,
  xlab = "1-SD tăng lên ở BMI so với \nmỗi 1-SD tăng lên ở nồng đồ của biomarker",
  title = "Mối liên quan giữa \ncác biomarkers trong máu và chỉ số BMI"
)