32 Đường cong dịch bệnh

Đường cong dịch bệnh (còn được gọi là “đường cong epi”) là một biểu đồ dịch tễ học cốt lõi thường được sử dụng để trực quan xu hướng khởi phát bệnh theo thời gian trong một cụm hoặc nhóm ca bệnh.

Phân tích đường cong dịch bệnh có thể cho biết xu hướng theo thời gian, giá trị ngoại lai, mức độ bùng phát dịch, khoảng thời gian có khả năng bị phơi nhiễm cao nhất, khoảng thời gian giữa các thế hệ ca bệnh và thậm chí có thể giúp xác định phương thức lây truyền của một căn bệnh không xác định (ví dụ: điểm bắt nguồn, nguồn tiếp diễn phổ biến, lây truyền từ người sang người). Bạn có thể tìm thấy một bài học trực tuyến về giải thích các đường cong dịch bệnh tại trang web của CDC Hoa Kỳ.

Trong chương này, chúng tôi trình bày hai cách tiếp cận để tạo ra đường cong dịch bệnh trong R:

  • Package incidence2, có thể tạo ra đường cong dịch bệnh với các lệnh đơn giản
  • Package ggplot2, cho phép khả năng tùy chỉnh nâng cao thông qua các lệnh phức tạp hơn

Chúng tôi cũng sẽ giải quyết các trường hợp cụ thể như:

  • Lập biểu đồ dữ liệu đếm tổng hợp
  • Faceting hoặc tạo nhiều cấu phần nhỏ
  • Áp dụng đường trung bình động
  • Hiển thị dữ liệu “dự kiến” hoặc có thể bị chậm trễ trong báo cáo
  • Thêm tỷ lệ ca nhiễm mới tích lũy bằng cách sử dụng trục thứ hai

32.1 Chuẩn bị

Packages

Đoạn code này hiển thị cách gọi các package cần thiết cho các phân tích. Trong cuốn sách này, chúng tôi nhấn mạnh đến hàm p_load() từ package pacman, sẽ cài đặt package nếu cần gọi package để sử dụng. Bạn cũng có thể cài đặt package với hàm library() từ base R. Xem chương R cơ bản để biết thêm thông tin về các package trong R.

pacman::p_load(
  rio,          # file import/export
  here,         # relative filepaths 
  lubridate,    # working with dates/epiweeks
  aweek,        # alternative package for working with dates/epiweeks
  incidence2,   # epicurves of linelist data
  i2extras,     # supplement to incidence2
  stringr,      # search and manipulate character strings
  forcats,      # working with factors
  RColorBrewer, # Color palettes from colorbrewer2.org
  tidyverse     # data management + ggplot2 graphics
) 

Nhập dữ liệu

Hai bộ dữ liệu mẫu được sử dụng trong chương này:

  • Bộ số liệu linelist về các ca bệnh từ một vụ dịch mô phỏng
  • Số lượng tổng hợp theo bệnh viện từ cùng một dịch bệnh mô phỏng bên trên

Các bộ dữ liệu được nhập bằng hàm import() từ package rio. Xem chương Nhập xuất dữ liệu để biết các cách nhập dữ liệu khác nhau.

Bộ số liệu linelist

Chúng ta nhập bộ dữ liệu về các ca bệnh mô phỏng từ một vụ dịch Ebola. Nếu bạn muốn tải xuống dữ liệu để làm theo từng bước, hãy xem hướng dẫn trong chương Tải sách và dữ liệu. Chúng tôi giả sử các tệp tin nằm trong thư mục làm việc nên không có thư mục con nào được chỉ định trong đường dẫn này.

linelist <- import("linelist_cleaned.xlsx")

50 dòng đầu tiên được hiển thị như bên dưới:

Số lượng ca bệnh do bệnh viện tổng hợp

Theo mục tiêu của cuốn sách này, bộ dữ liệu về số lượng tổng hợp hàng tuần theo bệnh viện được tạo từ linelist với đoạn code sau.

# import the counts data into R
count_data <- linelist %>% 
  group_by(hospital, date_hospitalisation) %>% 
  summarize(n_cases = dplyr::n()) %>% 
  filter(date_hospitalisation > as.Date("2013-06-01")) %>% 
  ungroup()

50 dòng đầu tiên được hiển thị như bên dưới:

Thiết lập các tham số

Để tạo báo cáo, bạn có thể muốn thiết lập các thông số có thể chỉnh sửa, chẳng hạn như ngày dữ liệu hiện tại (“ngày dữ liệu”). Sau đó, bạn có thể tham chiếu đối tượng data_date trong code khi áp dụng bộ lọc hoặc trong chú thích động.

## set the report date for the report
## note: can be set to Sys.Date() for the current date
data_date <- as.Date("2015-05-15")

Xác minh ngày

Xác minh rằng mỗi cột ngày có liên quan là phân lớp Ngày và có phạm vi giá trị thích hợp. Bạn có thể thực hiện việc này đơn giản bằng cách sử dụng hàm hist() cho histogram hoặc hàm range() với na.rm=TRUE, hoặc với hàm ggplot() như bên dưới.

# check range of onset dates
ggplot(data = linelist)+
  geom_histogram(aes(x = date_onset))

32.2 Đường cong dịch bệnh với package incidence2

Dưới đây, chúng tôi trình bày cách tạo các đường cong dịch bệnh bằng cách sử dụng package precision2. Các tác giả của package này đã cố gắng cho phép người dùng tạo và sửa đổi các đường cong dịch bệnh mà không cần biết cú pháp ggplot2. Phần lớn nội dung của chương này được điều chỉnh từ minh họa của package, bạn có thể tìm thấy tại trang github của package incidence2.

Ví dụ đơn giản

Cần có 2 bước để vẽ đường cong dịch bệnh với package incidence2:

  1. Tạo một incidence object (sử dụng hàm incidence())

    • Cung cấp dữ liệu
    • Xác định cột ngày tại date_index =
    • Xác định khoảng thời gian interval = các ca nên được tổng hợp (hàng ngày, tuần, tháng..)
    • Xác định bất kỳ cột dùng để nhóm nào (ví dụ như giới tính, bệnh viện, kết quả)
  2. Vẽ biểu đồ incidence object

    • Xác định nhãn, màu, tiêu đề, …

Dưới đây, chúng tôi gọi package incidence2, tạo incidence object từ dữ liệu linelist trên cột date_onset và tổng hợp các trường hợp theo ngày. Sau đó, chúng tôi in ra tóm tắt về đối tượng incidence object.

# load incidence2 package
pacman::p_load(incidence2)

# create the incidence object, aggregating cases by day
epi_day <- incidence(       # create incidence object
  x = linelist,             # dataset
  date_index = date_onset,  # date column
  interval = "day"          # date grouping interval
  )

Đối tượng incidence2 trông giống như một tibble (một kiểu data frame) và có thể được hiển thị hoặc thao tác thêm như một bộ dữ liệu.

class(epi_day)
## [1] "incidence2"   "incidence_df" "tbl_df"       "tbl"          "data.frame"

Đây là thông tin khi được hiển thị. Nó có một cột date_index và một cột count .

epi_day
## An incidence object: 367 x 2
## date range: [2014-04-07] to [2015-04-30]
## cases: 5632
## interval: 1 day
## 
##    date_index count
##    <date>     <int>
##  1 2014-04-07     1
##  2 2014-04-15     1
##  3 2014-04-21     2
##  4 2014-04-25     1
##  5 2014-04-26     1
##  6 2014-04-27     1
##  7 2014-05-01     2
##  8 2014-05-03     1
##  9 2014-05-04     1
## 10 2014-05-05     1
## # ... with 357 more rows

Bạn cũng có thể tổng hợp thông tin của đối tượng:

# print summary of the incidence object
summary(epi_day)
## date range: [2014-04-07] to [2015-04-30]
## cases: 5632
## interval: 1 day
## timespan: 389 days

Để vẽ biểu đồ incidence object, hãy sử dụng hàm plot() với tên của incidence object. Trong nền, hàm plot.incidence2() đã được gọi, vì vậy để đọc tài liệu cụ thể về incidence2, bạn cần chạy lệnh ?plot.incidence2.

# plot the incidence object
plot(epi_day)

Nếu bạn nhận thấy nhiều đường dọc nhỏ màu trắng, hãy cố gắng điều chỉnh kích thước hình ảnh của bạn. Ví dụ: nếu bạn xuất biểu đồ bằng ggsave(), bạn có thể cung cấp giá trị số cho width =height =. Nếu bạn mở rộng biểu đồ, những đường đó có thể biến mất.

Thay đổi khoảng thời gian tổng hợp ca bệnh

Đối số interval = của hàm incidence() xác định cách các quan sát được nhóm thành các cột dọc trong biểu đồ.

Xác định khoảng thời gian

incidence2 đem lại sự linh hoạt và cú pháp dễ hiểu để bạn có thể dễ dàng tổng hợp số ca bệnh thành các biểu đồ cột-đường cong dịch bệnh. Cung cấp một giá trị như những giá trị bên dưới cho đối số interval =. Bạn có thể viết bất kỳ thông tin nào bên dưới dưới dạng số nhiều (ví dụ: “weeks”) và bạn có thể thêm số phía trước (ví dụ: “3 months”).

Tùy chọn đối số Giải thích thêm
Số (1, 7, 13, 14, etc.) số ngày cho mỗi khoảng thời gian
“weeks” ghi chú: mặc định ngày bắt đầu là Thứ 2
“2 weeks” hoặc 3, 4, 5…
“Sunday weeks” tuần bắt đầu tính từ Chủ nhật (cũng có thể là Thứ 5, …)
“2 Sunday weeks” hoặc 3, 4, 5…
“MMWRweek” tuần bắt đầu vào Chủ nhật - xem thêm tại CDC Hoa Kỳ
“month” ngày đầu tiên của tháng
“quarter” ngày đầu tiên của tháng của một quí
“2 months” hoặc 3, 4, 5…
“year” ngày đầu tiên của năm

Dưới đây là ví dụ về các khoảng thời gian khác nhau trông như thế nào khi áp dụng vào bộ số liệu linelist. Lưu ý cách mà định dạng mặc định và tần suất của các nhãn thời gian trên trục x thay đổi khi khoảng thời gian thay đổi.

# Create the incidence objects (with different intervals)
##############################
# Weekly (Monday week by default)
epi_wk      <- incidence(linelist, date_onset, interval = "Monday week")

# Sunday week
epi_Sun_wk  <- incidence(linelist, date_onset, interval = "Sunday week")

# Three weeks (Monday weeks by default)
epi_2wk     <- incidence(linelist, date_onset, interval = "2 weeks")

# Monthly
epi_month   <- incidence(linelist, date_onset, interval = "month")

# Quarterly
epi_quarter   <- incidence(linelist, date_onset, interval = "quarter")

# Years
epi_year   <- incidence(linelist, date_onset, interval = "year")


# Plot the incidence objects (+ titles for clarity)
############################
plot(epi_wk)+      labs(title = "Monday weeks")
plot(epi_Sun_wk)+  labs(title = "Sunday weeks")
plot(epi_2wk)+     labs(title = "2 (Monday) weeks")
plot(epi_month)+   labs(title = "Months")
plot(epi_quarter)+ labs(title = "Quarters")
plot(epi_year)+    labs(title = "Years")

Ngày đầu tiên

Bạn có thể tùy chọn chỉ định một giá trị thuộc phân lớp Ngày (ví dụ: as.Date("2016-05-01")) tới đối số firstdate = trong lệnh incidence(). Nếu được cung cấp, dữ liệu sẽ được cắt bớt theo phạm vi này và khoảng thời gian sẽ bắt đầu từ ngày này.

Nhóm

Các nhóm được chỉ định trong lệnh incidence() và có thể được sử dụng để tô màu cho các cột hoặc facet dữ liệu. Để chỉ định các nhóm trong dữ liệu của bạn, hãy cung cấp tên cột tới đối số groups = trong hàm incidence() (không có dấu ngoặc kép xung quanh tên cột). Nếu chỉ định nhiều cột, hãy để tên các cột bên trong hàm c().

Bạn có thể chỉ định rằng các trường hợp có giá trị bị thiếu trong các cột được nhóm được liệt kê như một nhóm NA riêng biệt bằng cách thiết lập na_as_group = TRUE. Nếu không, chúng sẽ bị loại khỏi biểu đồ.

  • Để tô màu các cột theo cột nhóm, bạn cần cung cấp lại tên cột tới đối số fill = trong lệnh plot().

  • Để facet dựa trên cột nhóm, hãy xem phần bên dưới về facet với incidence2.

Trong ví dụ dưới đây, các trường hợp trong toàn bộ đợt bùng phát được nhóm theo phân loại tuổi của họ. Các giá trị thiếu được đưa vào một nhóm. Khoảng thời gian của đường cong bệnh dịch là tuần.

# Create incidence object, with data grouped by age category
age_outbreak <- incidence(
  linelist,                # dataset
  date_index = date_onset, # date column
  interval = "week",       # Monday weekly aggregation of cases
  groups = age_cat,        # age_cat is set as a group
  na_as_group = TRUE)      # missing values assigned their own group

# plot the grouped incidence object
plot(
  age_outbreak,             # incidence object with age_cat as group
  fill = age_cat)+          # age_cat is used for bar fill color (must have been set as a groups column above)
labs(fill = "Age Category") # change legend title from default "age_cat" (this is a ggplot2 modification)

MẸO: Thay đổi tiêu đề của chú giải bằng cách thêm + lệnh ggplot2 labs(fill = "your title") vào biểu đồ incidence2.

Bạn cũng có thể để các cột được nhóm hiển thị cạnh nhau bằng cách đặt stack = FALSE trong plot(), như được hiển thị bên dưới:

# Make incidence object of monthly counts. 
monthly_gender <- incidence(
 linelist,
 date_index = date_onset,
 interval = "month",
 groups = gender            # set gender as grouping column
)

plot(
  monthly_gender,   # incidence object
  fill = gender,    # display bars colored by gender
  stack = FALSE)    # side-by-side (not stacked)

Bạn có thể đặt đối số na_as_group = thành FALSE trong lệnh incidence() để loại bỏ các hàng có giá trị bị thiếu khỏi biểu đồ.

Lọc dữ liệu

Để vẽ biểu đồ đường cong dịch bệnh của một bộ dữ liệu con:

  1. Lọc bộ số liệu linelist
  2. Cung cấp dữ liệu đã lọc vào lệnh incidence()
  3. Vẽ biểu đồ incidence object

Ví dụ dưới đây sử dụng dữ liệu được lọc để chỉ hiển thị các trường hợp tại Central Hospital.

# filter the linelist
central_data <- linelist %>% 
  filter(hospital == "Central Hospital")

# create incidence object using filtered data
central_outbreak <- incidence(central_data, date_index = date_onset, interval = "week")

# plot the incidence object
plot(central_outbreak, title = "Weekly case incidence at Central Hospital")

Số lượng tổng hợp

Nếu dữ liệu ban đầu của bạn được tổng hợp (số lượng), hãy cung cấp tên của cột chứa thông tin về số lượng ca tới đối số count = khi tạo incidence object với hàm incidence().

Ví dụ, data frame count_data này được tổng hợp từ linelist theo số lượng hàng ngày theo bệnh viện. 50 hàng đầu tiên hiển thị như sau:

Nếu bạn đang bắt đầu phân tích với dữ liệu đếm hàng ngày như bộ dữ liệu ở trên, thì lệnh incidence() để chuyển đổi dữ liệu này thành đường cong dịch bệnh hàng tuần theo bệnh viện sẽ trông như sau:

epi_counts <- incidence(              # create weekly incidence object
  count_data,                         # dataset with counts aggregated by day
  date_index = date_hospitalisation,  # column with dates
  count = n_cases,                    # column with counts
  interval = "week",                  # aggregate daily counts up to weeks
  groups = hospital                   # group by hospital
  )

# plot the weekly incidence epi curve, with stacked bars by hospital
plot(epi_counts,                      # incidence object
     fill = hospital)                 # color the bars by hospital

Facets/Các biểu đồ nhỏ

Để facet dữ liệu theo nhóm (ví dụ tạo “các biểu đồ nhỏ”):

  1. Xác định các cột dùng để facet tới đối số groups = khi bạn tạo incidence object
  2. Sử dụng lệnh facet_plot() thay cho plot()
  3. Xác định cột nhóm nào được sử dụng cho fill = và cho facets =

Dưới đây, chúng ta thiết lập cả hai cột hospitaloutcome làm cột chia nhóm trong lệnh incidence(). Sau đó, trong lệnh facet_plot(), chúng ta vẽ biểu đồ đường cong dịch bệnh, chỉ rõ rằng chúng ta muốn có các đường cong dịch bệnh khác nhau cho mỗi bệnh viện và trong mỗi đường cong dịch bệnh, các cột được xếp chồng lên nhau và được tô màu theo outcome.

epi_wks_hosp_out <- incidence(
  linelist,                      # dataset
  date_index = date_onset,       # date column
  interval = "month",            # monthly bars  
  groups = c(outcome, hospital)  # both outcome and hospital are given as grouping columns
  )

# plot
incidence2::facet_plot(
  epi_wks_hosp_out,      # incidence object
  facets = hospital,     # facet column
  fill = outcome)        # fill column

Lưu ý rằng package ggtree (được sử dụng để trực quan cây phả hệ) cũng có một hàm là facet_plot() - đây là lý do tại sao chúng ta phải cụ thể incidence2::facet_plot() như ở trên.

Hiệu chỉnh với plot()

Một đường cong dịch bệnh được tạo ra bởi incidence2 có thể được sửa đổi thông qua các đối số trong hàm plot().

Dưới đây là các đối số của hàm plot() dùng để điều chỉnh hình thức cột:

Đối số Mô tả Ví dụ
fill = Màu cột. Tên màu hoặc tên cột đã được chỉ định trước đó tới đối số groups = trong lệnh incidence() fill = "red", hoặc fill = gender
color = Tô màu xung quanh mỗi cột hoặc xung quanh mỗi nhóm trong một cột border = "white"
legend = Vị trí danh mục Một trong các lựa chọn “bottom”, “top”, “left”, “right”, hoặc “none”
alpha = Độ trong suốt của cột 1 là đục hoàn toàn, 0 là trong suốt hoàn toàn
width = Giá trị từ 0 đến 1 cho biết kích thước tương đối của các cột với khoảng thời gian của chúng width = .7
show_cases = Giá trị logic; nếu TRUE, mỗi trường hợp hiển thị như một hộp. Cách hiển thị này là tốt nhất trên các đợt bùng phát nhỏ. show_cases = TRUE

Dưới đây là các đối số của hàm plot() để sửa đổi trục ngày tháng:

Đối số Mô tả
centre_dates = TRUE/FALSE tương đương với việc hiển thị ngày tháng xuất hiện ở dưới trung tâm các cột hay ở trên đỉnh các cột
date_format = Điều chỉnh định dạng hiển thị ngày bằng cú pháp strptime("%"). Chỉ hoạt động nếu centre_dates = FALSE (chi tiết bên dưới).
n.breaks = Khoảng số lượng nhãn ngắt khoảng trên trục x mong muốn.
angle = Góc nghiêng của nhãn ngày trên trục x (số độ)
size = Kích thước của văn bản tính bằng điểm

Lưu ý rằng đối số date_breaks = chỉ hoạt động nếu centre_dates = FALSE. Cung cấp giá trị ký tự trong dấu ngoặc kép bằng cú pháp strptime như bên dưới, như đã được trình bày chi tiết trong chương Làm việc với ngày tháng. Bạn có thể sử dụng \n để tạo “dòng mới”.

%d = Số ngày của tháng (5, 17, 28, …)
%j = Số ngày của năm (ngày Julian 001-366)
%a = Ngày trong tuần viết tắt (Mon, Tue, Wed, …)
%A = Ngày trong tuần viết đầy đủ (Monday, Tuesday, …)
%w = Số thứ tự ngày của tuần (0-6, Chủ nhật là 0)
%u = Số thứ tự ngày của tuần (1-7, Thứ hai là 1)
%W = Số thứ tự tuần (00-53, Thứ hai là ngày bắt đầu tuần mới)
%U = Số thứ tự tuần (01-53, Chủ nhật là ngày bắt đầu tuần mới)
%m = Số thứ tự tháng (Ví dụ 01, 02, 03, 04)
%b = Tháng viết tắt (Jan, Feb, …)
%B = Tháng viết đầy đủ (January, February, …)
%y = Năm viết dạng 2 ký tự (Ví dụ 89)
%Y = Năm viết dạng 2 ký tự (Ví dụ 1989)
%H = giờ (khung 24 giờ)
%M = phút
%S = giây
%z = bù từ múi giờ GMT
%Z = Múi giờ (ký tự)

Dưới đây là các đối số trong hàm plot() để hiệu chỉnh các nhãn trong biểu đồ:

Đối số Mô tả
title = Tiêu đề của biểu đồ
xlab = Tiêu đề của trục x
ylab = Tiêu đề của trục y
size = Kích thước của đoạn văn bản trên trục x tính bằng đơn vị pts (sử dụng hàm theme() trong ggplot để điều chỉnh các thông số kích thước khác)

Một ví dụ sử dụng kết hợp các đối số ở trên:

# filter the linelist
central_data <- linelist %>% 
  filter(hospital == "Central Hospital")

# create incidence object using filtered data
central_outbreak <- incidence(
  central_data,
  date_index = date_onset,
  interval = "week",
  groups = outcome)

# plot incidence object
plot(
  central_outbreak,
  fill = outcome,                       # box/bar color
  legend = "top",                       # legend on top
  title = "Cases at Central Hospital",  # title
  xlab = "Week of onset",               # x-axis label
  ylab = "Week of onset",               # y-axis label
  show_cases = TRUE,                    # show each case as an individual box
  alpha = 0.7,                          # transparency 
  border = "grey",                      # box border
  angle = 30,                           # angle of date labels
  centre_dates = FALSE,                 # date labels at edge of bar
  date_format = "%a %d %b %Y\n(Week %W)" # adjust how dates are displayed
  )

Để điều chỉnh thêm giao diện hiển thị của biểu đồ, hãy xem mục bên dưới về các hiệu chỉnh với hàm ggplot().

Hiệu chỉnh với ggplot2

Bạn có thể hiệu chỉnh biểu đồ được tạo bởi incidence2 bằng cách thêm các hiệu chỉnh ggplot2 với dấu + sau khi đóng ngoặc hàm biểu diễn tỷ lệ mắc bệnh plot(), như được minh họa bên dưới.

Dưới đây, biểu đồ incidence2 kết thúc và ngay sau đó các lệnh ggplot2 được sử dụng để sửa đổi các trục, thêm chú thích và điều chỉnh phông chữ đậm và kích thước chữ.

Lưu ý rằng nếu bạn thêm scale_x_date(), hầu hết các định dạng ngày từ plot() sẽ bị ghi đè. Xem mục đường cong dịch bệnh với ggplot() và chương Các tips với ggplot để có thêm lựa chọn.

# filter the linelist
central_data <- linelist %>% 
  filter(hospital == "Central Hospital")

# create incidence object using filtered data
central_outbreak <- incidence(
  central_data,
  date_index = date_onset,
  interval = "week",
  groups = c(outcome))

# plot incidence object
plot(
  central_outbreak,
  fill = outcome,                       # box/bar color
  legend = "top",                       # legend on top
  title = "Cases at Central Hospital",  # title
  xlab = "Week of onset",               # x-axis label
  ylab = "Week of onset",               # y-axis label
  show_cases = TRUE,                    # show each case as an individual box
  alpha = 0.7,                          # transparency 
  border = "grey",                      # box border
  centre_dates = FALSE,                   
  date_format = "%a %d %b\n%Y (Week %W)", 
  angle = 30                           # angle of date labels
  )+
  
  scale_y_continuous(
    breaks = seq(from = 0, to = 30, by = 5),  # specify y-axis increments by 5
    expand = c(0,0))+                         # remove excess space below 0 on y-axis
  
  # add dynamic caption
  labs(
    fill = "Patient outcome",                               # Legend title
    caption = stringr::str_glue(                            # dynamic caption - see page on characters and strings for details
      "n = {central_cases} from Central Hospital
      Case onsets range from {earliest_date} to {latest_date}. {missing_onset} cases are missing date of onset and not shown",
      central_cases = nrow(central_data),
      earliest_date = format(min(central_data$date_onset, na.rm=T), format = '%a %d %b %Y'),
      latest_date = format(max(central_data$date_onset, na.rm=T), format = '%a %d %b %Y'),      
      missing_onset = nrow(central_data %>% filter(is.na(date_onset)))))+
  
  # adjust bold face, and caption position
  theme(
    axis.title = element_text(size = 12, face = "bold"),    # axis titles larger and bold
    axis.text = element_text(size = 10, face = "bold"),     # axis text size and bold
    plot.caption = element_text(hjust = 0, face = "italic") # move caption to left
  )

Đổi màu

Chỉ định bảng màu

Cung cấp tên của bảng màu đã được định danh trước tới tham số col_pal = trong hàm plot(). Package incidence2 đi kèm với 2 bảng màu được xác định trước: “vibrant” và “muted”. Trong “vibrant”, có 6 màu đầu tiên là riêng biệt và trong “muted” có 9 màu đầu tiên là khác biệt. Sau những con số này, các màu là sự thêm vào/trung gian của các màu khác. Bạn có thể tìm thấy những bảng màu được xác định trước này tại website này. Các bảng màu loại trừ màu xám, do nó được dành riêng cho dữ liệu bị thiếu (sử dụng na_color = để thay đổi mặc định này).

# Create incidence object, with data grouped by age category  
age_outbreak <- incidence(
  linelist,
  date_index = date_onset,   # date of onset for x-axis
  interval = "week",         # weekly aggregation of cases
  groups = age_cat)

# plot the epicurve with default palette
plot(age_outbreak, fill = age_cat, title = "'vibrant' default incidence2 palette")

# plot with different color palette
#plot(age_outbreak, fill = age_cat, col_pal = muted, title = "'muted' incidence2 palette")

Bạn cũng có thể sử dụng một trong các bảng màu trong base R (đặt tên của bảng màu không có dấu ngoặc kép).

# plot with base R palette
plot(age_outbreak, fill = age_cat, col_pal = heat.colors, title = "base R heat.colors palette")

# plot with base R palette
plot(age_outbreak, fill = age_cat, col_pal = rainbow, title = "base R rainbow palette")

Bạn cũng có thể thêm bảng màu từ package viridis hoặc package RColorBrewer. Đầu tiên các package đó cần được tải, sau đó thêm các hàm scale_fill_*() tương ứng của chúng bằng dấu +, như bên dưới đây.

pacman::p_load(RColorBrewer, viridis)

# plot with color palette
plot(age_outbreak, fill = age_cat, title = "Viridis palette")+
  scale_fill_viridis_d(
    option = "inferno",     # color scheme, try also "plasma" or the default
    name = "Age Category",  # legend name
    na.value = "grey")      # for missing values

# plot with color palette
plot(age_outbreak, fill = age_cat, title = "RColorBrewer palette")+
  scale_fill_brewer(
    palette = "Dark2",      # color palette, try also Accent, Dark2, Paired, Pastel1, Pastel2, Set1, Set2, Set3
    name = "Age Category",  # legend name
    na.value = "grey")      # for missing values

Xác định thủ công

Để chỉ định màu theo cách thủ công, hãy thêm hàm scale_fill_manual() trong ggplot2 sau hàm plot() bằng dấu + và cung cấp vectơ tên màu hoặc mã HEX tới đối số values =. Số lượng màu được liệt kê phải bằng số nhóm. Hãy lưu ý, các giá trị bị thiếu có phải là một nhóm hay không - chúng có thể được chuyển đổi thành giá trị dạng ký tự chẳng hạn như “Missing” trong quá trình chuẩn bị dữ liệu của bạn bằng hàm fct_explicit_na(), như đã được giải thích trong chương Factors.

# manual colors
plot(age_outbreak, fill = age_cat, title = "Manually-specified colors")+
  scale_fill_manual(
    values = c("darkgreen", "darkblue", "purple", "grey", "yellow", "orange", "red", "lightblue"),  # colors
    name = "Age Category")      # Name for legend

Như đã đề cập trong chương Các tips với ggplot, bạn có thể tạo bảng màu của riêng mình bằng cách sử dụng hàm colorRampPalette() trên một vectơ chứa tên các màu sắc và chỉ định số màu bạn muốn trả về. Đây là một cách tốt để có thể thu được nhiều màu trong một lần bằng cách chỉ định một vài màu.

my_cols <- c("darkgreen", "darkblue", "purple", "grey", "yellow", "orange")
my_palette <- colorRampPalette(my_cols)(12)  # expand the 6 colors above to 12 colors
my_palette
##  [1] "#006400" "#00363F" "#00097E" "#3A0BAF" "#821ADD" "#A84BE2" "#B592CB" "#C9C99B" "#E7E745" "#FFF600" "#FFCD00" "#FFA500"

Điều chỉnh thứ bậc

Để điều chỉnh thứ tự xuất hiện của nhóm (trên biểu đồ và trong chú giải), cột phân nhóm cần phải là kiểu Factor. Xem chương Factors để biết thêm thông tin.

Đầu tiên, chúng ta hãy quan sát một đường cong dịch bệnh hàng tuần theo bệnh viện với thứ tự mặc định:

# ORIGINAL - hospital NOT as factor
###################################

# create weekly incidence object, rows grouped by hospital and week
hospital_outbreak <- incidence(
  linelist,
  date_index = date_onset, 
  interval = "week", 
  groups = hospital)

# plot incidence object
plot(hospital_outbreak, fill = hospital, title = "ORIGINAL - hospital not a factor")

Bây giờ, để điều chỉnh thứ tự sao cho nhóm phân loại “Missing” và “Other” ở trên cùng của đường cong dịch bệnh, chúng ta có thể làm như sau:

  • Gọi package forcats, để làm việc với cột định dạng factor

  • Điều chỉnh bộ dữ liệu - trong trường hợp này, chúng tôi sẽ xác định một bộ dữ liệu mới (plot_data) trong đó:

    • Cột gender được xác định là một biến phân loại có thứ bậc được thiết lập bằng hàm fct_relevel() sao cho “Other” và “Missing” đứng đầu tiên, do đó chúng xuất hiện ở trên đầu các cột
  • Incidence object được tạo và vẽ biểu đồ như trước đó

  • Chúng ta thêm các điều chỉnh ggplot2

    • scale_fill_manual() để gán màu theo cách thủ công sao cho “Missing” là màu xám và “Other” là màu be
# MODIFIED - hospital as factor
###############################

# load forcats package for working with factors
pacman::p_load(forcats)

# Convert hospital column to factor and adjust levels
plot_data <- linelist %>% 
  mutate(hospital = fct_relevel(hospital, c("Missing", "Other"))) # Set "Missing" and "Other" as top levels


# Create weekly incidence object, grouped by hospital and week
hospital_outbreak_mod <- incidence(
  plot_data,
  date_index = date_onset, 
  interval = "week", 
  groups = hospital)

# plot incidence object
plot(hospital_outbreak_mod, fill = hospital)+
  
  # manual specify colors
  scale_fill_manual(values = c("grey", "beige", "darkgreen", "green2", "orange", "red", "pink"))+                      

  # labels added via ggplot
  labs(
      title = "MODIFIED - hospital as factor",   # plot title
      subtitle = "Other & Missing at top of epicurve",
      y = "Weekly case incidence",               # y axis title  
      x = "Week of symptom onset",               # x axis title
      fill = "Hospital")                         # title of legend     

MẸO: Nếu bạn chỉ muốn đảo ngược thứ tự của chú giải, hãy thêm lệnh ggplot2 này guides(fill = guide_legend(reverse = TRUE)).

Đường lưới dọc

Nếu bạn vẽ biểu đồ với các thiết lập mặc định của incidence2, bạn có thể nhận thấy rằng các đường lưới dọc xuất hiện ở từng nhãn ngày và giữa mỗi nhãn ngày. Điều này có thể dẫn đến các đường lưới giao nhau với đỉnh của một số cột.

Bạn có thể xóa tất cả các đường lưới bằng cách thêm lệnh ggplot2 theme_classic().

# make incidence object
a <- incidence(
  central_data,
  date_index = date_onset,
  interval = "Monday weeks"
)

# Default gridlines
plot(a, title = "Default lines")

# Specified gridline intervals
# NOT WORKING WITH INCIDENCE2 1.0.0
# plot(a, title = "Weekly lines")+
#   scale_x_date(
#     date_breaks = "4 weeks",      # major vertical lines align on weeks
#     date_minor_breaks = "weeks",  # minor vertical lines every week
#     date_labels = "%a\n%d\n%b")   # format of date labels

# No gridlines
plot(a, title = "No lines")+
  theme_classic()                 # remove all gridlines

Tuy nhiên, lưu ý rằng nếu sử dụng tuần, các đối số date_breaksdate_minor_breaks chỉ hoạt động cho các tuần được quy định bắt đầu vào Thứ Hai. Nếu các tuần của bạn bắt đầu vào một ngày khác trong tuần, thay vào đó, bạn sẽ cần cung cấp vectơ ngày cho các đối số breaks =minor_breaks = theo cách thủ công. Xem mục ggplot2 để biết các ví dụ về điều này bằng cách sử dụng hàm seq.Date().

Số mới mắc tích lũy

Xem phần sau để biết phương pháp thay thế để vẽ biểu đồ số mới mắc tích lũy với ggplot2 - ví dụ: để chồng một đường số mới mắc tích lũy lên trên đường cong dịch bệnh.

Trung bình động

Bạn có thể thêm đường trung bình động vào biểu đồ tạo bởi incidence2 một cách dễ dàng với hàm add_rolling_average() từ package i2extras. Chuyển đối tượng incidence2 của bạn tới hàm này, và sau đó vẽ bằng hàm plot(). Thiết lập before = là số ngày trước đó bạn muốn đưa vào đường trung bình động (mặc định là 2). Nếu dữ liệu của bạn được nhóm lại, trung bình động sẽ được tính cho mỗi nhóm.

rolling_avg <- incidence(                    # make incidence object
  linelist,
  date_index = date_onset,
  interval = "week",
  groups = gender) %>% 
  
  i2extras::add_rolling_average(before = 6)  # add rolling averages (in this case, by gender)

# plot
plot(rolling_avg, n.breaks = 3) # faceted automatically because rolling average on groups

Để tìm hiểu cách áp dụng đường trung bình động một cách tổng quát hơn trên dữ liệu, hãy xem chương Đường trung bình động.

32.3 Đường cong dịch bệnh với ggplot2

Sử dụng ggplot() để xây dựng đường cong dịch bệnh cho phép sự linh hoạt và tùy chỉnh nhiều hơn, nhưng đòi hỏi nhiều nỗ lực và hiểu biết hơn về sử dụng ggplot().

Không giống như sử dụng package incidence2, bạn phải kiểm soát thủ công việc tổng hợp các trường hợp theo thời gian (theo tuần, tháng, v.v.) khoảng thời gian của các nhãn trên trục ngày thág. Việc này phải được quản lý rất cẩn thận.

Các ví dụ này sử dụng một tập con của bộ dữ liệu linelist - chỉ các trường hợp từ Central Hospital.

central_data <- linelist %>% 
  filter(hospital == "Central Hospital")

Để tạo ra một đường cong dịch bệnh với ggplot(), có ba yếu tố chính:

  • Một histogram, với các ca bệnh linelist được tổng hợp thành các “bins” được phân biệt bằng các điểm “ngắt” cụ thể
  • Cân chỉnh các trục và nhãn
  • Chủ đề của biểu đồ, bao gồm tiêu đề, nhãn, chú thích, v.v.

Xác định bins của số liệu

Ở đây chúng ta trình bày cách xác định các ca bệnh sẽ được tổng hợp thành các bins của một histogram (“cột”). Điều quan trọng cần nhận ra là việc tổng hợp số ca bệnh thành các bins của một histogram không nhất thiết phải có cùng khoảng với thời gian xuất hiện trên trục x.

Dưới đây có lẽ là đoạn code đơn giản nhất để tạo ra các đường cong dịch bệnh hàng ngày và hàng tuần.

Trong lệnh ggplot() tổng quát, bộ dữ liệu được đưa vào data =. Trên nền tảng này, dạng hình học của histogram được thêm bằng một dấu +. Trong geom_histogram(), chúng ta ánh xạ các yếu tố thẩm mỹ sao cho cột date_onset được ánh xạ tới trục x. Cũng trong hàm geom_histogram() nhưng không phải trong đối số aes(), chúng ta đặt binwidth = bằng số bins của histogram, tính bằng ngày. Nếu cú pháp ggplot2 này khó hiểu, hãy xem lại chương ggplot cơ bản.

THẬN TRỌNG: Lập biểu đồ các trường hợp hàng tuần bằng cách sử dụng binwidth = 7 bắt đầu cột 7 ngày đầu tiên ở ca bệnh đầu tiên, mà có thể là bất kỳ ngày nào trong tuần! Để tạo các tuần cụ thể, hãy xem phần bên dưới.

# daily 
ggplot(data = central_data) +          # set data
  geom_histogram(                      # add histogram
    mapping = aes(x = date_onset),     # map date column to x-axis
    binwidth = 1)+                     # cases binned by 1 day 
  labs(title = "Central Hospital - Daily")                # title

# weekly
ggplot(data = central_data) +          # set data 
  geom_histogram(                      # add histogram
      mapping = aes(x = date_onset),   # map date column to x-axis
      binwidth = 7)+                   # cases binned every 7 days, starting from first case (!) 
  labs(title = "Central Hospital - 7-day bins, starting at first case") # title

Chúng tôi xin lưu ý rằng ca bệnh đầu tiên trong bộ dữ liệu của Central Hospital có triệu chứng khởi phát vào ngày:

format(min(central_data$date_onset, na.rm=T), "%A %d %b, %Y")
## [1] "Thursday 01 May, 2014"

Để chỉ định thủ công chia biểu đồ cột, không sử dụng argument binwidth = và thay vào đó cung cấp vectơ ngày để breaks =.

Tạo vectơ ngày tháng với hàm seq.Date() trong base R. Hàm này sử dụng các đối số to =, from =, và by =. Ví dụ: lệnh bên dưới trả về ngày mỗi tháng bắt đầu từ ngày 15 Tháng 1 và kết thúc trước ngày 28 Tháng 6.

monthly_breaks <- seq.Date(from = as.Date("2014-02-01"),
                           to = as.Date("2015-07-15"),
                           by = "months")

monthly_breaks   # print
##  [1] "2014-02-01" "2014-03-01" "2014-04-01" "2014-05-01" "2014-06-01" "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01" "2014-12-01"
## [12] "2015-01-01" "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01" "2015-07-01"

Vectơ này có được cung cấp cho hàm geom_histogram() dưới dạng breaks =:

# monthly 
ggplot(data = central_data) +  
  geom_histogram(
    mapping = aes(x = date_onset),
    breaks = monthly_breaks)+         # provide the pre-defined vector of breaks                    
  labs(title = "Monthly case bins")   # title

Một chuỗi ngày hàng tuần đơn giản có thể được trả về bằng cách đặt by = "week". Ví dụ:

weekly_breaks <- seq.Date(from = as.Date("2014-02-01"),
                          to = as.Date("2015-07-15"),
                          by = "week")

Một giải pháp thay thế cho việc xác định ngày bắt đầu và ngày kết thúc cụ thể là viết code động để các cột hàng tuần bắt đầu vào thứ Hai trước ca đầu tiên. Chúng tôi sẽ sử dụng các vectơ ngày này trong suốt các ví dụ dưới đây.

# Sequence of weekly Monday dates for CENTRAL HOSPITAL
weekly_breaks_central <- seq.Date(
  from = floor_date(min(central_data$date_onset, na.rm=T),   "week", week_start = 1), # monday before first case
  to   = ceiling_date(max(central_data$date_onset, na.rm=T), "week", week_start = 1), # monday after last case
  by   = "week")

Hãy giải mã đoạn code khá “khoai” ở trên:

  • Giá trị “from” (ngày sớm nhất của chuỗi) được tạo như sau: giá trị ngày nhỏ nhất (min() với na.rm=TRUE) trong cột date_onset được đưa vào hàm floor_date() thuộc package lubridate. floor_date() được đặt thành “tuần” trả về ngày bắt đầu của “tuần” của trường hợp đó, với điều kiện là ngày bắt đầu của mỗi tuần là Thứ Hai (week_start = 1).
  • Tương tự như vậy, giá trị “to” (ngày kết thúc của chuỗi) được tạo bằng cách sử dụng hàm ngược lại ceiling_date() để trả về thứ hai sau ca cuối cùng.
  • Đối số “by” của seq.Date() có thể được đặt thành bất kỳ số ngày, tuần hoặc tháng nào.
  • Sử dụng week_start = 7 cho tuần bắt đầu vào Chủ nhật

Vì chúng ta sẽ sử dụng các vectơ ngày dạng này trong toàn bộ chương này, chúng ta cũng xác định một vectơ cho toàn bộ đợt bùng phát (ở trên chỉ dành cho Central Hospital).

# Sequence for the entire outbreak
weekly_breaks_all <- seq.Date(
  from = floor_date(min(linelist$date_onset, na.rm=T),   "week", week_start = 1), # monday before first case
  to   = ceiling_date(max(linelist$date_onset, na.rm=T), "week", week_start = 1), # monday after last case
  by   = "week")

Các kết quả trả về của seq.Date() có thể được sử dụng để tạo các khoảng chia bins trong histogram, cũng như khoảng chia cho các nhãn ngày mà có thể độc lập với các bins. Đọc thêm về nhãn ngày trong các phần sau.

MẸO: Đối với lệnh ggplot() đơn giản hơn, hãy lưu các khoảng chia bins và khoảng chia nhãn ngày dưới dạng các vectơ đã đặt tên trước và gán tên này vào đối số breaks =.

Ví dụ về đường cong dịch bệnh theo tuần

Dưới đây là code ví dụ chi tiết để tạo đường cong dịch bệnh hàng tuần cho các tuần bắt đầu vào Thứ Hai, với các cột, nhãn ngày và đường lưới dọc đã được căn chỉnh. Phần này dành cho người dùng có nhu cầu code nhanh. Để hiểu sâu từng khía cạnh (chủ đề, nhãn ngày, v.v.), hãy tiếp tục xem các phần tiếp theo. Chú ý:

  • Các đoạn chia bins của histogram được xác định bằng hàm seq.Date() như được giải thích ở trên để bắt đầu vào thứ Hai trước ca sớm nhất và kết thúc vào thứ Hai sau ca cuối cùng
  • Khoảng nhãn ngày được xác định bởi date_breaks = bên trong scale_x_date()
  • Khoảng đường lưới dọc nhỏ giữa các nhãn ngày được xác định bởi date_minor_breaks =
  • expand = c(0,0) trong trục x và y loại bỏ không gian thừa trên mỗi cạnh của trục, điều này cũng đảm bảo các nhãn ngày bắt đầu từ cột đầu tiên.
# TOTAL MONDAY WEEK ALIGNMENT
#############################
# Define sequence of weekly breaks
weekly_breaks_central <- seq.Date(
      from = floor_date(min(central_data$date_onset, na.rm=T),   "week", week_start = 1), # Monday before first case
      to   = ceiling_date(max(central_data$date_onset, na.rm=T), "week", week_start = 1), # Monday after last case
      by   = "week")    # bins are 7-days 


ggplot(data = central_data) + 
  
  # make histogram: specify bin break points: starts the Monday before first case, end Monday after last case
  geom_histogram(
    
    # mapping aesthetics
    mapping = aes(x = date_onset),  # date column mapped to x-axis
    
    # histogram bin breaks
    breaks = weekly_breaks_central, # histogram bin breaks defined previously
    
    # bars
    color = "darkblue",     # color of lines around bars
    fill = "lightblue"      # color of fill within bars
  )+ 
    
  # x-axis labels
  scale_x_date(
    expand            = c(0,0),           # remove excess x-axis space before and after case bars
    date_breaks       = "4 weeks",        # date labels and major vertical gridlines appear every 3 Monday weeks
    date_minor_breaks = "week",           # minor vertical lines appear every Monday week
    date_labels       = "%a\n%d %b\n%Y")+ # date labels format
  
  # y-axis
  scale_y_continuous(
    expand = c(0,0))+             # remove excess y-axis space below 0 (align histogram flush with x-axis)
  
  # aesthetic themes
  theme_minimal()+                # simplify plot background
  
  theme(
    plot.caption = element_text(hjust = 0,        # caption on left side
                                face = "italic"), # caption in italics
    axis.title = element_text(face = "bold"))+    # axis titles in bold
  
  # labels including dynamic caption
  labs(
    title    = "Weekly incidence of cases (Monday weeks)",
    subtitle = "Note alignment of bars, vertical gridlines, and axis labels on Monday weeks",
    x        = "Week of symptom onset",
    y        = "Weekly incident cases reported",
    caption  = stringr::str_glue("n = {nrow(central_data)} from Central Hospital; Case onsets range from {format(min(central_data$date_onset, na.rm=T), format = '%a %d %b %Y')} to {format(max(central_data$date_onset, na.rm=T), format = '%a %d %b %Y')}\n{nrow(central_data %>% filter(is.na(date_onset)))} cases missing date of onset and not shown"))

Tuần bắt đầu bằng Chủ nhật

Để vẽ được biểu đồ cho các tuần bắt đầu bằng Chủ nhật, cần có một số sửa đổi, vì date_breaks = "weeks" chỉ hoạt động cho các tuần bắt đầu bằng thứ Hai.

  • Các điểm ngắt của các histogram bins phải được đặt thành Chủ nhật (week_start = 7)
  • Trong scale_x_date(), các đoạn ngắt ngày tương tự nên được gắn với breaks =minor_breaks = để đảm bảo các nhãn ngày và đường lưới dọc căn chỉnh vào ngày Chủ nhật.

Ví dụ: lệnh scale_x_date() cho các tuần bắt đầu vào Chủ nhật có thể trông như sau:

scale_x_date(
    expand = c(0,0),
    
    # specify interval of date labels and major vertical gridlines
    breaks = seq.Date(
      from = floor_date(min(central_data$date_onset, na.rm=T),   "week", week_start = 7), # Sunday before first case
      to   = ceiling_date(max(central_data$date_onset, na.rm=T), "week", week_start = 7), # Sunday after last case
      by   = "4 weeks"),
    
    # specify interval of minor vertical gridline 
    minor_breaks = seq.Date(
      from = floor_date(min(central_data$date_onset, na.rm=T),   "week", week_start = 7), # Sunday before first case
      to   = ceiling_date(max(central_data$date_onset, na.rm=T), "week", week_start = 7), # Sunday after last case
      by   = "week"),
   
    # date label format
    date_labels = "%a\n%d %b\n%Y")+         # day, above month abbrev., above 2-digit year

Nhóm/tô màu theo giá trị

Các cột histogram có thể được tô màu theo nhóm và “xếp chồng lên nhau”. Để chỉ định cột nhóm, hãy thực hiện các thay đổi sau. Xem chương ggplot cơ bản để biết thêm chi tiết.

  • Bên trong đối số aes() của hàm vẽ histogram, ánh xạ tên cột tới đối số group =fill =
  • Loại bỏ bất kỳ tham số fill = nào bên ngoài đối số aes(), vì nó sẽ ghi đè giá trị bên trong
  • Đối số bên trong aes() sẽ áp dụng theo nhóm, trong khi bất kỳ đối số bên ngoài nào đều sẽ áp dụng cho tất cả các cột (ví dụ: bạn có thể vẫn muốn tô màu color = bên ngoài, vì vậy mỗi cột có cùng một đường viền)

Đây là lệnh aes() hiển thị nhóm và tô màu các cột theo giới tính:

aes(x = date_onset, group = gender, fill = gender)

Đây là kết quả khi áp dụng:

ggplot(data = linelist) +     # begin with linelist (many hospitals)
  
  # make histogram: specify bin break points: starts the Monday before first case, end Monday after last case
  geom_histogram(
    mapping = aes(
      x = date_onset,
      group = hospital,       # set data to be grouped by hospital
      fill = hospital),       # bar fill (inside color) by hospital
    
    # bin breaks are Monday weeks
    breaks = weekly_breaks_all,   # sequence of weekly Monday bin breaks for whole outbreak, defined in previous code       
    
    # Color around bars
    color = "black")

Hiệu chỉnh màu

  • Để thiết lập thủ công tô màu cho từng nhóm, hãy sử dụng scale_fill_manual() (lưu ý: scale_color_manual() là một cái khác!)

    • Sử dụng đối số values = để áp dụng một vectơ màu
    • Sử dụng na.value = để xác định màu cho giá trị NA
    • Sử dụng đối số labels = để thay đổi văn bản trong mục chú thích. Để cho an toàn, sử dụng một vectơ được đặt tên kiểu như c("old" = "new", "old" = "new") hoặc điều chỉnh giá trị trong bộ dữ liệu
    • Sử dụng name = để đặt một tiêu đề thích hợp cho mục chú thích
  • Để biết thêm mẹo về thang màu và bảng màu, hãy xem chương ggplot cơ bản.

ggplot(data = linelist)+           # begin with linelist (many hospitals)
  
  # make histogram
  geom_histogram(
    mapping = aes(x = date_onset,
        group = hospital,          # cases grouped by hospital
        fill = hospital),          # bar fill by hospital
    
    # bin breaks
    breaks = weekly_breaks_all,        # sequence of weekly Monday bin breaks, defined in previous code
    
    # Color around bars
    color = "black")+              # border color of each bar
  
  # manual specification of colors
  scale_fill_manual(
    values = c("black", "orange", "grey", "beige", "blue", "brown"),
    labels = c("St. Mark's Maternity Hospital (SMMH)" = "St. Mark's"),
    name = "Hospital") # specify fill colors ("values") - attention to order!

Hiệu chỉnh thứ bậc

Thứ tự mà các cột được nhóm xếp chồng lên nhau được điều chỉnh tốt nhất bằng cách phân loại cột nhóm dưới dạng Factor. Sau đó, bạn có thể chỉ định thứ tự cấp độ phân loại (và nhãn hiển thị của chúng). Xem chương về Factors hoặc Các tips với ggplot để biết thêm chi tiết.

Trước khi tạo biểu đồ, hãy sử dụng hàm fct_relevel() từ package forcats để chuyển đổi cột phân nhóm thành kiểu Factor và điều chỉnh thứ tự cấp độ theo cách thủ công, như được trình bày chi tiết trong chương về Factors.

# load forcats package for working with factors
pacman::p_load(forcats)

# Define new dataset with hospital as factor
plot_data <- linelist %>% 
  mutate(hospital = fct_relevel(hospital, c("Missing", "Other"))) # Convert to factor and set "Missing" and "Other" as top levels to appear on epicurve top

levels(plot_data$hospital) # print levels in order
## [1] "Missing"                              "Other"                                "Central Hospital"                    
## [4] "Military Hospital"                    "Port Hospital"                        "St. Mark's Maternity Hospital (SMMH)"

Trong biểu đồ dưới đây, điểm khác biệt duy nhất so với trước đó là cột hospital đã được hợp nhất như trên và chúng ta sử dụng các guides() để đảo ngược thứ tự chú thích, do đó “nhãn Missing” nằm ở cuối chú thích.

ggplot(plot_data) +                     # Use NEW dataset with hospital as re-ordered factor
  
  # make histogram
  geom_histogram(
    mapping = aes(x = date_onset,
        group = hospital,               # cases grouped by hospital
        fill = hospital),               # bar fill (color) by hospital
    
    breaks = weekly_breaks_all,         # sequence of weekly Monday bin breaks for whole outbreak, defined at top of ggplot section
    
    color = "black")+                   # border color around each bar
    
  # x-axis labels
  scale_x_date(
    expand            = c(0,0),         # remove excess x-axis space before and after case bars
    date_breaks       = "3 weeks",      # labels appear every 3 Monday weeks
    date_minor_breaks = "week",         # vertical lines appear every Monday week
    date_labels       = "%d\n%b\n'%y")+ # date labels format
  
  # y-axis
  scale_y_continuous(
    expand = c(0,0))+                   # remove excess y-axis space below 0
  
  # manual specification of colors, ! attention to order
  scale_fill_manual(
    values = c("grey", "beige", "black", "orange", "blue", "brown"),
    labels = c("St. Mark's Maternity Hospital (SMMH)" = "St. Mark's"),
    name = "Hospital")+ 
  
  # aesthetic themes
  theme_minimal()+                      # simplify plot background
  
  theme(
    plot.caption = element_text(face = "italic", # caption on left side in italics
                                hjust = 0), 
    axis.title = element_text(face = "bold"))+   # axis titles in bold
  
  # labels
  labs(
    title    = "Weekly incidence of cases by hospital",
    subtitle = "Hospital as re-ordered factor",
    x        = "Week of symptom onset",
    y        = "Weekly cases")

MẸO: Để chỉ đảo ngược thứ tự của chú thích, hãy thêm lệnh ggplot2 này: guides(fill = guide_legend(reverse = TRUE)).

Hiệu chỉnh chú thích

Đọc thêm về chú thích và scale trong chương Các tips với ggplot. Dưới đây là một vài điểm nổi bật:

  • Chỉnh sửa tiêu đề chú thích trong hàm scale hoặc với labs(fill = "Legend title") (nếu bạn đang sử dụng color =, thì hãy sử dụng labs(color = ""))

  • theme(legend.title = element_blank()) để bỏ trống tiêu đề chú thích

  • theme(legend.position = "top") (“bottom”, “left”, “right”, hoặc “none” để bỏ chú thích)

  • theme(legend.direction = "horizontal") xoay ngang chú thích

  • guides(fill = guide_legend(reverse = TRUE)) để đảo ngược thứ tự các mục chú thích

Cột kề cột

Hiển thị song song các cột nhóm (trái ngược với xếp chồng lên nhau) được chỉ định trong geom_histogram() với position = "dodge" được đặt bên ngoài aes().

Nếu có nhiều hơn hai nhóm giá trị, nó có thể gây khó đọc. Thay vào đó, hãy cân nhắc sử dụng một biểu đồ được chia nhỏ (gồm nhiều biểu đồ nhỏ). Để dễ xem trong ví dụ này, các giá trị giới tính bị thiếu sẽ bị xóa.

ggplot(central_data %>% drop_na(gender))+   # begin with Central Hospital cases dropping missing gender
    geom_histogram(
        mapping = aes(
          x = date_onset,
          group = gender,         # cases grouped by gender
          fill = gender),         # bars filled by gender
        
        # histogram bin breaks
        breaks = weekly_breaks_central,   # sequence of weekly dates for Central outbreak - defined at top of ggplot section
        
        color = "black",          # bar edge color
        
        position = "dodge")+      # SIDE-BY-SIDE bars
                      
  
  # The labels on the x-axis
  scale_x_date(expand            = c(0,0),         # remove excess x-axis space below and after case bars
               date_breaks       = "3 weeks",      # labels appear every 3 Monday weeks
               date_minor_breaks = "week",         # vertical lines appear every Monday week
               date_labels       = "%d\n%b\n'%y")+ # date labels format
  
  # y-axis
  scale_y_continuous(expand = c(0,0))+             # removes excess y-axis space between bottom of bars and the labels
  
  #scale of colors and legend labels
  scale_fill_manual(values = c("brown", "orange"),  # specify fill colors ("values") - attention to order!
                    na.value = "grey" )+     

  # aesthetic themes
  theme_minimal()+                                               # a set of themes to simplify plot
  theme(plot.caption = element_text(face = "italic", hjust = 0), # caption on left side in italics
        axis.title = element_text(face = "bold"))+               # axis titles in bold
  
  # labels
  labs(title    = "Weekly incidence of cases, by gender",
       subtitle = "Subtitle",
       fill     = "Gender",                                      # provide new title for legend
       x        = "Week of symptom onset",
       y        = "Weekly incident cases reported")

Giới hạn trục

Có hai cách để giới hạn phạm vi các giá trị trên trục.

Nói chung, cách khuyến khích sử dụng là lệnh coord_cartesian(), chấp nhận xlim = c(min, max)ylim = c(min, max) (trong đó bạn cung cấp giá trị nhỏ nhất và lớn nhất). Cách này hoạt động như một “thu phóng kích thước” mà không thực sự loại bỏ bất kỳ dữ liệu nào, điều này rất quan trọng đối với các số liệu thống kê và các thang đo tổng hợp.

Ngoài ra, bạn có thể thiết lập giá trị ngày tối đa và tối thiểu bằng cách sử dụng limits = c() bên trong hàm scale_x_date(). Ví dụ:

scale_x_date(limits = c(as.Date("2014-04-01"), NA)) # sets a minimum date but leaves the maximum open.  

Tương tự như vậy, nếu bạn muốn trục x kéo dài đến một ngày cụ thể (ví dụ: ngày hiện tại), ngay cả khi không có trường hợp mới nào được báo cáo, bạn có thể sử dụng:

scale_x_date(limits = c(NA, Sys.Date()) # ensures date axis will extend until current date  

LƯU Ý: Hãy thận trọng khi đặt giới hạn hoặc chia tỷ lệ trục y (ví dụ: 0 đến 30 với khoảng là 5: seq(0, 30, 5)). Những con số tĩnh như vậy có thể cắt bớt biểu đồ quá ngắn nếu dữ liệu thay đổi vượt quá giới hạn!.

Nhãn/đường lưới của trục ngày

MẸO: Hãy nhớ rằng các nhãn trục ngày độc lập với việc tổng hợp dữ liệu thành các cột, nhưng về mặt trực quan, điều quan trọng là phải căn chỉnh các cột, nhãn ngày và các đường lưới dọc.

Để sửa đổi nhãn ngày và đường lưới, hãy sử dụng scale_x_date() theo một trong những cách sau:

  • Nếu histogram bins là ngày, tuần bắt đầu vào thứ hai, tháng hoặc năm:

    • Sử dụng date_breaks = để xác định khoảng thời gian của các nhãn và đường lưới chính (ví dụ: “day”, “week”, “3 weeks”, “month” hoặc “year”)
    • Sử dụng date_minor_breaks = để xác định khoảng của các đường lưới dọc nhỏ (giữa các nhãn ngày)
    • Thêm expand = c(0,0) để bắt đầu các nhãn ở cột đầu tiên
    • Sử dụng date_labels = để xác định định dạng của nhãn ngày - hãy xem chương Ngày để biết các mẹo (sử dụng \n cho một dòng mới)
  • Nếu histogram bins là các tuần bắt đầu vào Chủ nhật:

    • Sử dụng breaks =minor_breaks = bằng cách cung cấp một chuỗi ngày cho mỗi khoảng chia
    • Vẫn có thể sử dụng date_labels =expand = để định dạng như mô tả ở trên

Một số ghi chú:

Minh họa

Dưới đây là minh họa biểu đồ trong đó các cột và nhãn biểu đồ/đường lưới được căn chỉnh thẳng hàng và không thẳng hàng:

# 7-day bins + Monday labels
#############################
ggplot(central_data) +
  geom_histogram(
    mapping = aes(x = date_onset),
    binwidth = 7,                 # 7-day bins with start at first case
    color = "darkblue",
    fill = "lightblue") +
  
  scale_x_date(
    expand = c(0,0),               # remove excess x-axis space below and after case bars
    date_breaks = "3 weeks",       # Monday every 3 weeks
    date_minor_breaks = "week",    # Monday weeks
    date_labels = "%a\n%d\n%b\n'%y")+  # label format
  
  scale_y_continuous(
    expand = c(0,0))+              # remove excess space under x-axis, make flush
  
  labs(
    title = "MISALIGNED",
    subtitle = "! CAUTION: 7-day bars start Thursdays at first case\nDate labels and gridlines on Mondays\nNote how ticks don't align with bars")



# 7-day bins + Months
#####################
ggplot(central_data) +
  geom_histogram(
    mapping = aes(x = date_onset),
    binwidth = 7,
    color = "darkblue",
    fill = "lightblue") +
  
  scale_x_date(
    expand = c(0,0),                  # remove excess x-axis space below and after case bars
    date_breaks = "months",           # 1st of month
    date_minor_breaks = "week",       # Monday weeks
    date_labels = "%a\n%d %b\n%Y")+    # label format
  
  scale_y_continuous(
    expand = c(0,0))+                # remove excess space under x-axis, make flush 
  
  labs(
    title = "MISALIGNED",
    subtitle = "! CAUTION: 7-day bars start Thursdays with first case\nMajor gridlines and date labels at 1st of each month\nMinor gridlines weekly on Mondays\nNote uneven spacing of some gridlines and ticks unaligned with bars")


# TOTAL MONDAY ALIGNMENT: specify manual bin breaks to be mondays
#################################################################
ggplot(central_data) + 
  geom_histogram(
    mapping = aes(x = date_onset),
    
    # histogram breaks set to 7 days beginning Monday before first case
    breaks = weekly_breaks_central,    # defined earlier in this page
    
    color = "darkblue",
    
    fill = "lightblue") + 
  
  scale_x_date(
    expand = c(0,0),                   # remove excess x-axis space below and after case bars
    date_breaks = "4 weeks",           # Monday every 4 weeks
    date_minor_breaks = "week",        # Monday weeks 
    date_labels = "%a\n%d %b\n%Y")+      # label format
  
  scale_y_continuous(
    expand = c(0,0))+                # remove excess space under x-axis, make flush 
  
  labs(
    title = "ALIGNED Mondays",
    subtitle = "7-day bins manually set to begin Monday before first case (28 Apr)\nDate labels and gridlines on Mondays as well")


# TOTAL MONDAY ALIGNMENT WITH MONTHS LABELS:
############################################
ggplot(central_data) + 
  geom_histogram(
    mapping = aes(x = date_onset),
    
    # histogram breaks set to 7 days beginning Monday before first case
    breaks = weekly_breaks_central,            # defined earlier in this page
    
    color = "darkblue",
    
    fill = "lightblue") + 
  
  scale_x_date(
    expand = c(0,0),                   # remove excess x-axis space below and after case bars
    date_breaks = "months",            # Monday every 4 weeks
    date_minor_breaks = "week",        # Monday weeks 
    date_labels = "%b\n%Y")+          # label format
  
  scale_y_continuous(
    expand = c(0,0))+                # remove excess space under x-axis, make flush 
  
  theme(panel.grid.major = element_blank())+  # Remove major gridlines (fall on 1st of month)
          
  labs(
    title = "ALIGNED Mondays with MONTHLY labels",
    subtitle = "7-day bins manually set to begin Monday before first case (28 Apr)\nDate labels on 1st of Month\nMonthly major gridlines removed")


# TOTAL SUNDAY ALIGNMENT: specify manual bin breaks AND labels to be Sundays
############################################################################
ggplot(central_data) + 
  geom_histogram(
    mapping = aes(x = date_onset),
    
    # histogram breaks set to 7 days beginning Sunday before first case
    breaks = seq.Date(from = floor_date(min(central_data$date_onset, na.rm=T),   "week", week_start = 7),
                      to   = ceiling_date(max(central_data$date_onset, na.rm=T), "week", week_start = 7),
                      by   = "7 days"),
    
    color = "darkblue",
    
    fill = "lightblue") + 
  
  scale_x_date(
    expand = c(0,0),
    # date label breaks and major gridlines set to every 3 weeks beginning Sunday before first case
    breaks = seq.Date(from = floor_date(min(central_data$date_onset, na.rm=T),   "week", week_start = 7),
                      to   = ceiling_date(max(central_data$date_onset, na.rm=T), "week", week_start = 7),
                      by   = "3 weeks"),
    
    # minor gridlines set to weekly beginning Sunday before first case
    minor_breaks = seq.Date(from = floor_date(min(central_data$date_onset, na.rm=T),   "week", week_start = 7),
                            to   = ceiling_date(max(central_data$date_onset, na.rm=T), "week", week_start = 7),
                            by   = "7 days"),
    
    date_labels = "%a\n%d\n%b\n'%y")+  # label format
  
  scale_y_continuous(
    expand = c(0,0))+                # remove excess space under x-axis, make flush 
  
  labs(title = "ALIGNED Sundays",
       subtitle = "7-day bins manually set to begin Sunday before first case (27 Apr)\nDate labels and gridlines manually set to Sundays as well")

Dữ liệu tổng hợp

Thông thường, thay vì bắt đầu với bộ số liệu linelist, bạn có thể bắt đầu với số lượng tổng hợp từ các cơ sở y tế, quận, huyện, v.v. Bạn có thể tạo đường cong dịch bệnh với ggplot() nhưng code sẽ hơi khác một chút. Phần này sẽ tận dụng bộ dữ liệu count_data đã được nạp trước đó trong mục chuẩn bị dữ liệu. Bộ dữ liệu này là linelist được tổng hợp thành số lượng bệnh viện theo ngày. 50 hàng đầu tiên được hiển thị dưới đây.

Vẽ biểu đồ dữ liệu đếm hàng ngày

Chúng ta có thể vẽ biểu đồ đường cong dịch bệnh từ dữ liệu đếm hàng ngày. Dưới đây là sự khác biệt trong đoạn code:

  • Khi ánh xạ các yếu tố trong hàm aes(), hãy cụ thể y = bằng một cột đếm (trong trường hợp này, tên cột là n_cases)

  • Thêm đối số stat = "identity" trong geom_histogram(), giúp xác định rằng chiều cao cột phải là giá trị y =, không phải số hàng như mặc định

  • Thêm đối số width = để tránh các đường trắng dọc giữa các cột. Đối với dữ liệu hàng ngày được đặt thành 1. Đối với dữ liệu đếm hàng tuần được đặt thành 7. Đối với dữ liệu đếm hàng tháng, các đường màu trắng là một vấn đề (mỗi tháng có số ngày khác nhau) - hãy xem xét chuyển đổi trục x của bạn thành cột phân loại (factor) theo thứ tự (tháng) và sử dụng geom_col().

ggplot(data = count_data)+
  geom_histogram(
   mapping = aes(x = date_hospitalisation, y = n_cases),
   stat = "identity",
   width = 1)+                # for daily counts, set width = 1 to avoid white space between bars
  labs(
    x = "Date of report", 
    y = "Number of cases",
    title = "Daily case incidence, from daily count data")

Vẽ biểu đồ dữ liệu đếm hàng tuần

Nếu dữ liệu của bạn đã là số lượng ca bệnh đếm theo tuần, chúng có thể trông giống như bộ dữ liệu này (được gọi là count_data_weekly):

50 hàng đầu tiên của count_data_weekly được hiển thị bên dưới. Bạn có thể thấy rằng số lượng đếm đã được tổng hợp thành các tuần. Mỗi tuần được hiển thị theo ngày đầu tiên của tuần (mặc định là Thứ Hai).

Bây giờ vẽ biểu đồ sao cho x = cột epiweek. Hãy nhớ thêm y = cột đếm khi ánh xạ trục và thêm stat = "identity" như đã giải thích ở trên.

ggplot(data = count_data_weekly)+
  
  geom_histogram(
    mapping = aes(
      x = epiweek,           # x-axis is epiweek (as class Date)
      y = n_cases_weekly,    # y-axis height in the weekly case counts
      group = hospital,      # we are grouping the bars and coloring by hospital
      fill = hospital),
    stat = "identity")+      # this is also required when plotting count data
     
  # labels for x-axis
  scale_x_date(
    date_breaks = "2 months",      # labels every 2 months 
    date_minor_breaks = "1 month", # gridlines every month
    date_labels = '%b\n%Y')+       #labeled by month with year below
     
  # Choose color palette (uses RColorBrewer package)
  scale_fill_brewer(palette = "Pastel2")+ 
  
  theme_minimal()+
  
  labs(
    x = "Week of onset", 
    y = "Weekly case incidence",
    fill = "Hospital",
    title = "Weekly case incidence, from aggregated count data by hospital")

Đường trung bình động

Xem chương về Đường trung bình động để có mô tả chi tiết và một số tùy chọn. Dưới đây là một tùy chọn để tính toán đường trung bình động với package slider. Theo cách tiếp cận này, trung bình động được tính toán trong bộ dữ liệu trước khi vẽ biểu đồ:

  1. Tổng hợp dữ liệu thành số lượng nếu cần thiết (hàng ngày, hàng tuần, v.v.) (xem chương Nhóm dữ liệu)
  2. Tạo một cột mới để giữ đường trung bình động, được tạo bằng hàm slide_index() từ package slider
  3. Vẽ đường trung bình động dưới dạng một geom_line() ở trên đỉnh (phía sau) histogram đường cong dịch bệnh

Tham khảo thêm tại vignette của package slider

# load package
pacman::p_load(slider)  # slider used to calculate rolling averages

# make dataset of daily counts and 7-day moving average
#######################################################
ll_counts_7day <- linelist %>%    # begin with linelist
  
  ## count cases by date
  count(date_onset, name = "new_cases") %>%   # name new column with counts as "new_cases"
  drop_na(date_onset) %>%                     # remove cases with missing date_onset
  
  ## calculate the average number of cases in 7-day window
  mutate(
    avg_7day = slider::slide_index(    # create new column
      new_cases,                       # calculate based on value in new_cases column
      .i = date_onset,                 # index is date_onset col, so non-present dates are included in window 
      .f = ~mean(.x, na.rm = TRUE),    # function is mean() with missing values removed
      .before = 6,                     # window is the day and 6-days before
      .complete = FALSE),              # must be FALSE for unlist() to work in next step
    avg_7day = unlist(avg_7day))       # convert class list to class numeric


# plot
######
ggplot(data = ll_counts_7day) +  # begin with new dataset defined above 
    geom_histogram(              # create epicurve histogram
      mapping = aes(
        x = date_onset,          # date column as x-axis
        y = new_cases),          # height is number of daily new cases
        stat = "identity",       # height is y value
        fill="#92a8d1",          # cool color for bars
        colour = "#92a8d1",      # same color for bar border
        )+ 
    geom_line(                   # make line for rolling average
      mapping = aes(
        x = date_onset,          # date column for x-axis
        y = avg_7day,            # y-value set to rolling average column
        lty = "7-day \nrolling avg"), # name of line in legend
      color="red",               # color of line
      size = 1) +                # width of line
    scale_x_date(                # date scale
      date_breaks = "1 month",
      date_labels = '%d/%m',
      expand = c(0,0)) +
    scale_y_continuous(          # y-axis scale
      expand = c(0,0),
      limits = c(0, NA)) +       
    labs(
      x="",
      y ="Number of confirmed cases",
      fill = "Legend")+ 
    theme_minimal()+
    theme(legend.title = element_blank())  # removes title of legend

Faceting/chia nhỏ biểu đồ

Như với các ggplots khác, bạn có thể tạo các biểu đồ được chia nhỏ (“nhiều biểu đồ con”). Như đã giải thích trong chương Các tips với ggplot trong cuốn sách này, bạn có thể sử dụng facet_wrap() hoặc facet_grid(). Ở đây chúng tôi sẽ minh họa bằng hàm facet_wrap(). Đối với đường cong dịch bệnh, sử dụng facet_wrap() thường dễ dàng hơn vì khả năng bạn thường chỉ cần chia nhỏ biểu đồ theo một biến.

Cú pháp chung là facet_wrap(rows ~ cols), trong đó bên trái dấu ngã (~) là tên của biến sẽ được trải trên các “hàng” của biểu đồ chia nhỏ và ở bên phải dấu ngã là tên của biến sẽ được trải trên các “cột” của biểu đồ chia nhỏ. Đơn giản nhất, chỉ cần sử dụng một tên cột, ở bên phải dấu ngã: facet_wrap(~age_cat).

Trục tự do
Bạn sẽ cần phải quyết định xem tỷ lệ của các trục cho mỗi biểu đồ nhỏ là “cố định” với cùng một kích thước (mặc định) hay “tự do” (nghĩa là chúng sẽ thay đổi dựa trên dữ liệu của chúng). Thực hiện điều này với đối số scales = trong hàm facet_wrap() bằng cách chỉ định “free_x” hoặc “free_y” hoặc “free”.

Số cột và hàng của các biểu đồ con
Điều này có thể xác định với ncol =nrow = trong hàm facet_wrap().

Thứ tự các biểu đồ con
Để thay đổi thứ tự xuất hiện, hãy thay đổi thứ tự cơ bản của các cấp của cột phân loại được sử dụng để tạo các biểu đồ con.

Đinh dạng trục
Kích thước phông chữ và mặt biểu đồ, dải màu, v.v. có thể được sửa đổi thông qua theme() với các đối số như:

  • strip.text = element_text() (kích thước, màu sắc, mặt, góc …)

  • strip.background = element_rect() (ví dụ: element_rect(fill = "grey"))

  • strip.position = (vị trị “dưới”, “trên”, “trái”, hoặc “phải”)

Dải nhãn
Nhãn của các biểu đồ con có thể được sửa đổi thông qua “nhãn” của cột như một factor hoặc bằng cách sử dụng một “người dán nhãn - labeller”.

Để tạo một labeller như thế, sử dụng hàm as_labeller() từ ggplot2. Sau đó, cung cấp labeller cho đối số labeller = của facet_wrap() như dưới đây.

my_labels <- as_labeller(c(
     "0-4"   = "Ages 0-4",
     "5-9"   = "Ages 5-9",
     "10-14" = "Ages 10-14",
     "15-19" = "Ages 15-19",
     "20-29" = "Ages 20-29",
     "30-49" = "Ages 30-49",
     "50-69" = "Ages 50-69",
     "70+"   = "Over age 70"))

Một ví dụ về chia nhỏ biểu đồ - chia bằng cột age_cat.

# make plot
###########
ggplot(central_data) + 
  
  geom_histogram(
    mapping = aes(
      x = date_onset,
      group = age_cat,
      fill = age_cat),    # arguments inside aes() apply by group
      
    color = "black",      # arguments outside aes() apply to all data
        
    # histogram breaks
    breaks = weekly_breaks_central)+  # pre-defined date vector (see earlier in this page)
                      
  # The labels on the x-axis
  scale_x_date(
    expand            = c(0,0),         # remove excess x-axis space below and after case bars
    date_breaks       = "2 months",     # labels appear every 2 months
    date_minor_breaks = "1 month",      # vertical lines appear every 1 month 
    date_labels       = "%b\n'%y")+     # date labels format
  
  # y-axis
  scale_y_continuous(expand = c(0,0))+                       # removes excess y-axis space between bottom of bars and the labels
  
  # aesthetic themes
  theme_minimal()+                                           # a set of themes to simplify plot
  theme(
    plot.caption = element_text(face = "italic", hjust = 0), # caption on left side in italics
    axis.title = element_text(face = "bold"),
    legend.position = "bottom",
    strip.text = element_text(face = "bold", size = 10),
    strip.background = element_rect(fill = "grey"))+         # axis titles in bold
  
  # create facets
  facet_wrap(
    ~age_cat,
    ncol = 4,
    strip.position = "top",
    labeller = my_labels)+             
  
  # labels
  labs(
    title    = "Weekly incidence of cases, by age category",
    subtitle = "Subtitle",
    fill     = "Age category",                                      # provide new title for legend
    x        = "Week of symptom onset",
    y        = "Weekly incident cases reported",
    caption  = stringr::str_glue("n = {nrow(central_data)} from Central Hospital; Case onsets range from {format(min(central_data$date_onset, na.rm=T), format = '%a %d %b %Y')} to {format(max(central_data$date_onset, na.rm=T), format = '%a %d %b %Y')}\n{nrow(central_data %>% filter(is.na(date_onset)))} cases missing date of onset and not shown"))

Xem link để biết thêm thông tin về labellers.

Tổng vụ dịch trong nền của biểu đồ con

Để hiển thị tổng vụ dịch trong nền của mỗi biểu đồ con, hãy thêm hàm gghighlight() với dấu ngoặc đơn trống vào ggplot. Hàm này thuộc package gghighlight. Lưu ý rằng trục y tối đa trong tất cả các biểu đồ con hiện dựa vào đỉnh của toàn bộ vụ dịch. Có nhiều ví dụ hơn về package này trong chương Các tips với ggplot.

ggplot(central_data) + 
  
  # epicurves by group
  geom_histogram(
    mapping = aes(
      x = date_onset,
      group = age_cat,
      fill = age_cat),  # arguments inside aes() apply by group
    
    color = "black",    # arguments outside aes() apply to all data
    
    # histogram breaks
    breaks = weekly_breaks_central)+     # pre-defined date vector (see top of ggplot section)                
  
  # add grey epidemic in background to each facet
  gghighlight::gghighlight()+
  
  # labels on x-axis
  scale_x_date(
    expand            = c(0,0),         # remove excess x-axis space below and after case bars
    date_breaks       = "2 months",     # labels appear every 2 months
    date_minor_breaks = "1 month",      # vertical lines appear every 1 month 
    date_labels       = "%b\n'%y")+     # date labels format
  
  # y-axis
  scale_y_continuous(expand = c(0,0))+  # removes excess y-axis space below 0
  
  # aesthetic themes
  theme_minimal()+                                           # a set of themes to simplify plot
  theme(
    plot.caption = element_text(face = "italic", hjust = 0), # caption on left side in italics
    axis.title = element_text(face = "bold"),
    legend.position = "bottom",
    strip.text = element_text(face = "bold", size = 10),
    strip.background = element_rect(fill = "white"))+        # axis titles in bold
  
  # create facets
  facet_wrap(
    ~age_cat,                          # each plot is one value of age_cat
    ncol = 4,                          # number of columns
    strip.position = "top",            # position of the facet title/strip
    labeller = my_labels)+             # labeller defines above
  
  # labels
  labs(
    title    = "Weekly incidence of cases, by age category",
    subtitle = "Subtitle",
    fill     = "Age category",                                      # provide new title for legend
    x        = "Week of symptom onset",
    y        = "Weekly incident cases reported",
    caption  = stringr::str_glue("n = {nrow(central_data)} from Central Hospital; Case onsets range from {format(min(central_data$date_onset, na.rm=T), format = '%a %d %b %Y')} to {format(max(central_data$date_onset, na.rm=T), format = '%a %d %b %Y')}\n{nrow(central_data %>% filter(is.na(date_onset)))} cases missing date of onset and not shown"))

Dữ liệu một biểu đồ con

Nếu bạn muốn có một hộp biểu đồ con chứa tất cả dữ liệu, hãy sao chép toàn bộ dữ liệu và coi các bản sao như một giá trị cho các biểu đồ con. Hàm “trợ giúp” CreateAllFacet() bên dưới có thể hỗ trợ việc này (nhờ bài đăng trên blog này). Khi nó được chạy, số hàng tăng gấp đôi và sẽ có một cột mới được gọi là facet, trong đó các hàng được sao chép sẽ có giá trị “tất cả” và các hàng ban đầu có giá trị ban đầu của cột phân chia. Bây giờ bạn chỉ cần phân chia cột facet .

Dưới đây là hàm trợ giúp. Chạy code này và nó sẽ luôn sẵn sàng để bạn sử dụng.

# Define helper function
CreateAllFacet <- function(df, col){
     df$facet <- df[[col]]
     temp <- df
     temp$facet <- "all"
     merged <-rbind(temp, df)
     
     # ensure the facet value is a factor
     merged[[col]] <- as.factor(merged[[col]])
     
     return(merged)
}

Bây giờ hãy áp dụng hàm trợ giúp cho bộ dữ liệu, trên cột age_cat:

# Create dataset that is duplicated and with new column "facet" to show "all" age categories as another facet level
central_data2 <- CreateAllFacet(central_data, col = "age_cat") %>%
  
  # set factor levels
  mutate(facet = fct_relevel(facet, "all", "0-4", "5-9",
                             "10-14", "15-19", "20-29",
                             "30-49", "50-69", "70+"))
## Warning: Unknown levels in `f`: 70+
# check levels
table(central_data2$facet, useNA = "always")
## 
##   all   0-4   5-9 10-14 15-19 20-29 30-49 50-69  <NA> 
##   454    84    84    82    58    73    57     7     9

Những thay đổi đáng chú ý đối với lệnh ggplot() là:

  • Dữ liệu được sử dụng bây giờ là central_data2 (nhân đôi các hàng, với cột mới là “facet”)

  • Labeller sẽ cần được cập nhật, nếu được sử dụng

  • Tùy chọn: để có các biểu đồ con xếp chồng lên nhau theo chiều dọc: cột chia được chuyển sang các hàng bên cạnh của phương trình và ở bên phải được thay thế bằng “.” (facet_wrap(facet~.)), và ncol = 1. Bạn cũng có thể cần điều chỉnh chiều rộng và chiều cao của ảnh biểu đồ đã lưu dưới dạng png (xem ggsave() trong chương Các tips với ggplot).

ggplot(central_data2) + 
  
  # actual epicurves by group
  geom_histogram(
        mapping = aes(
          x = date_onset,
          group = age_cat,
          fill = age_cat),  # arguments inside aes() apply by group
        color = "black",    # arguments outside aes() apply to all data
        
        # histogram breaks
        breaks = weekly_breaks_central)+    # pre-defined date vector (see top of ggplot section)
                     
  # Labels on x-axis
  scale_x_date(
    expand            = c(0,0),         # remove excess x-axis space below and after case bars
    date_breaks       = "2 months",     # labels appear every 2 months
    date_minor_breaks = "1 month",      # vertical lines appear every 1 month 
    date_labels       = "%b\n'%y")+     # date labels format
  
  # y-axis
  scale_y_continuous(expand = c(0,0))+  # removes excess y-axis space between bottom of bars and the labels
  
  # aesthetic themes
  theme_minimal()+                                           # a set of themes to simplify plot
  theme(
    plot.caption = element_text(face = "italic", hjust = 0), # caption on left side in italics
    axis.title = element_text(face = "bold"),
    legend.position = "bottom")+               
  
  # create facets
  facet_wrap(facet~. ,                            # each plot is one value of facet
             ncol = 1)+            

  # labels
  labs(title    = "Weekly incidence of cases, by age category",
       subtitle = "Subtitle",
       fill     = "Age category",                                      # provide new title for legend
       x        = "Week of symptom onset",
       y        = "Weekly incident cases reported",
       caption  = stringr::str_glue("n = {nrow(central_data)} from Central Hospital; Case onsets range from {format(min(central_data$date_onset, na.rm=T), format = '%a %d %b %Y')} to {format(max(central_data$date_onset, na.rm=T), format = '%a %d %b %Y')}\n{nrow(central_data %>% filter(is.na(date_onset)))} cases missing date of onset and not shown"))

32.4 Dữ liệu dự kiến

Dữ liệu gần đây nhất được biểu thị trong đường cong dịch tễ nên được đánh dấu là dự kiến hoặc có thể báo cáo chậm trễ. Điều này có thể được thực hiện bằng cách thêm một đường thẳng đứng và/hoặc hình chữ nhật trong một số ngày cụ thể. Đây là hai tùy chọn:

  1. Sử dụng annotate():

    • Để sử dụng dạng đường annotate(geom = "segment"). Cung cấp x, xend, y, và yend. Hiệu chỉnh kích thước, kiểu dòng (lty), và màu.
    • Để sử dụng dạng hình chữ nhật annotate(geom = "rect"). Cung cấp xmin/xmax/ymin/ymax. Hiệu chỉnh màu và hệ số alpha.
  2. Nhóm dữ liệu theo trạng thái tạm thời và tô màu các cột đó theo cách khác nhau

THẬN TRỌNG: Bạn có thể thử hàm geom_rect() để vẽ hình chữ nhật, nhưng việc điều chỉnh độ trong suốt không khả thi trong bối cảnh cảnh bộ số liệu linelist. Hàm này sẽ phủ lên một hình chữ nhật cho mỗi hàng/quan sát!. Sử dụng hệ số alpha rất thấp (ví dụ: 0.01) hoặc một cách tiếp cận khác.

Sử dụng annotate()

  • Trong annotate(geom = "rect"), đối số xminxmax cần được định dạng phân lớp ngày
  • Lưu ý rằng vì những dữ liệu này được tổng hợp thành các cột hàng tuần và cột cuối cùng kéo dài đến Thứ Hai sau điểm dữ liệu cuối cùng, vùng được tô bóng có thể bao gồm 4 tuần
  • Đây là một ví dụ trực tuyến về annotate()
ggplot(central_data) + 
  
  # histogram
  geom_histogram(
    mapping = aes(x = date_onset),
    
    breaks = weekly_breaks_central,   # pre-defined date vector - see top of ggplot section
    
    color = "darkblue",
    
    fill = "lightblue") +

  # scales
  scale_y_continuous(expand = c(0,0))+
  scale_x_date(
    expand = c(0,0),                   # remove excess x-axis space below and after case bars
    date_breaks = "1 month",           # 1st of month
    date_minor_breaks = "1 month",     # 1st of month
    date_labels = "%b\n'%y")+          # label format
  
  # labels and theme
  labs(
    title = "Using annotate()\nRectangle and line showing that data from last 21-days are tentative",
    x = "Week of symptom onset",
    y = "Weekly case indicence")+ 
  theme_minimal()+
  
  # add semi-transparent red rectangle to tentative data
  annotate(
    "rect",
    xmin  = as.Date(max(central_data$date_onset, na.rm = T) - 21), # note must be wrapped in as.Date()
    xmax  = as.Date(Inf),                                          # note must be wrapped in as.Date()
    ymin  = 0,
    ymax  = Inf,
    alpha = 0.2,          # alpha easy and intuitive to adjust using annotate()
    fill  = "red")+
  
  # add black vertical line on top of other layers
  annotate(
    "segment",
    x     = max(central_data$date_onset, na.rm = T) - 21, # 21 days before last data
    xend  = max(central_data$date_onset, na.rm = T) - 21, 
    y     = 0,         # line begins at y = 0
    yend  = Inf,       # line to top of plot
    size  = 2,         # line size
    color = "black",
    lty   = "solid")+   # linetype e.g. "solid", "dashed"

  # add text in rectangle
  annotate(
    "text",
    x = max(central_data$date_onset, na.rm = T) - 15,
    y = 15,
    label = "Subject to reporting delays",
    angle = 90)

Đường thẳng đứng màu đen có thể tạo ra với code bên dưới, nhưng sử dụng geom_vline(), bạn sẽ mất khả năng kiểm soát chiều cao:

geom_vline(xintercept = max(central_data$date_onset, na.rm = T) - 21,
           size = 2,
           color = "black")

Màu cột

Một cách tiếp cận thay thế có thể là điều chỉnh màu sắc hoặc cách hiển thị của chính các cột dữ liệu dự kiến. Bạn có thể tạo một cột mới trong giai đoạn chuẩn bị dữ liệu và sử dụng cột đó để nhóm dữ liệu, sao cho aes(fill = ) của dữ liệu tạm thời có thể có màu hoặc hệ số alpha khác với các cột khác.

# add column
############
plot_data <- central_data %>% 
  mutate(tentative = case_when(
    date_onset >= max(date_onset, na.rm=T) - 7 ~ "Tentative", # tenative if in last 7 days
    TRUE                                       ~ "Reliable")) # all else reliable

# plot
######
ggplot(plot_data, aes(x = date_onset, fill = tentative)) + 
  
  # histogram
  geom_histogram(
    breaks = weekly_breaks_central,   # pre-defined data vector, see top of ggplot page
    color = "black") +

  # scales
  scale_y_continuous(expand = c(0,0))+
  scale_fill_manual(values = c("lightblue", "grey"))+
  scale_x_date(
    expand = c(0,0),                   # remove excess x-axis space below and after case bars
    date_breaks = "3 weeks",           # Monday every 3 weeks
    date_minor_breaks = "week",        # Monday weeks 
    date_labels = "%d\n%b\n'%y")+      # label format
  
  # labels and theme
  labs(title = "Show days that are tentative reporting",
    subtitle = "")+ 
  theme_minimal()+
  theme(legend.title = element_blank())                 # remove title of legend

32.5 Nhãn ngày nhiều cấp độ

Nếu bạn muốn phân nhiều cấp nhãn ngày (ví dụ: tháng và năm) mà không sao chép các cấp nhãn cấp độ thấp hơn, hãy xem xét một trong các cách tiếp cận bên dưới:

Hãy nhớ rằng - bạn có thể sử dụng các công cụ như \n trong các đối số date_labels hoặc labels để đặt các phần của mỗi nhãn trên một dòng mới bên dưới. Tuy nhiên, đoạn code dưới đây giúp bạn thực hiện nhiều năm hoặc tháng (ví dụ) ở dòng thấp hơn và chỉ một lần. Một số lưu ý về code bên dưới:

  • Số lượng ca bệnh được tổng hợp thành các tuần vì lý do thẩm mỹ. Xem chương Đường cong dịch bệnh (tab dữ liệu tổng hợp) để biết chi tiết.
  • Một miền geom_area() được sử dụng thay vì một histogram, vì phương pháp tiếp cận chia biểu đồ dưới đây không hoạt động tốt với histogram.

Tổng hợp dữ đếm hàng tuần

# Create dataset of case counts by week
#######################################
central_weekly <- linelist %>%
  filter(hospital == "Central Hospital") %>%   # filter linelist
  mutate(week = lubridate::floor_date(date_onset, unit = "weeks")) %>%  
  count(week) %>%                              # summarize weekly case counts
  drop_na(week) %>%                            # remove cases with missing onset_date
  complete(                                    # fill-in all weeks with no cases reported
    week = seq.Date(
      from = min(week),   
      to   = max(week),
      by   = "week"),
    fill = list(n = 0))                        # convert new NA values to 0 counts

Vẽ biểu đồ

# plot with box border on year
##############################
ggplot(central_weekly) +
  geom_area(aes(x = week, y = n),    # make line, specify x and y
            stat = "identity") +             # because line height is count number
  scale_x_date(date_labels="%b",             # date label format show month 
               date_breaks="month",          # date labels on 1st of each month
               expand=c(0,0)) +              # remove excess space on each end
  scale_y_continuous(
    expand  = c(0,0))+                       # remove excess space below x-axis
  facet_grid(~lubridate::year(week), # facet on year (of Date class column)
             space="free_x",                
             scales="free_x",                # x-axes adapt to data range (not "fixed")
             switch="x") +                   # facet labels (year) on bottom
  theme_bw() +
  theme(strip.placement = "outside",         # facet labels placement
        strip.background = element_rect(fill = NA, # facet labels no fill grey border
                                        colour = "grey50"),
        panel.spacing = unit(0, "cm"))+      # no space between facet panels
  labs(title = "Nested year labels, grey label border")
# plot with no box border on year
#################################
ggplot(central_weekly,
       aes(x = week, y = n)) +              # establish x and y for entire plot
  geom_line(stat = "identity",              # make line, line height is count number
            color = "#69b3a2") +            # line color
  geom_point(size=1, color="#69b3a2") +     # make points at the weekly data points
  geom_area(fill = "#69b3a2",               # fill area below line
            alpha = 0.4)+                   # fill transparency
  scale_x_date(date_labels="%b",            # date label format show month 
               date_breaks="month",         # date labels on 1st of each month
               expand=c(0,0)) +             # remove excess space
  scale_y_continuous(
    expand  = c(0,0))+                      # remove excess space below x-axis
  facet_grid(~lubridate::year(week),        # facet on year (of Date class column)
             space="free_x",                
             scales="free_x",               # x-axes adapt to data range (not "fixed")
             switch="x") +                  # facet labels (year) on bottom
  theme_bw() +
  theme(strip.placement = "outside",                     # facet label placement
          strip.background = element_blank(),            # no facet lable background
          panel.grid.minor.x = element_blank(),          
          panel.border = element_rect(colour="grey40"),  # grey border to facet PANEL
          panel.spacing=unit(0,"cm"))+                   # No space between facet panels
  labs(title = "Nested year labels - points, shaded, no label border")

Các kỹ thuật trên được điều chỉnh từ đây và bài đăng này trên stackoverflow.com.

32.6 Trục kép

Mặc dù có những cuộc thảo luận gay gắt về tính hợp lệ của trục kép trong cộng đồng về trực quan hóa dữ liệu, nhiều chuyên gia dịch tễ vẫn muốn nhìn nhận biểu đồ đường cong dịch bệnh hoặc biểu đồ tương tự với phần trăm trên trục thứ hai. Điều này được thảo luận nhiều hơn trong chương Các tips với ggplot, nhưng một ví dụ sử dụng phương pháp cowplot được trình bày bên dưới:

  • Hai biểu đồ riêng biệt được tạo, và sau đó được kết hợp với package cowplot.
  • Các biểu đồ phải có chính xác cùng trục x (đã đặt giới hạn), nếu không dữ liệu và nhãn sẽ không được căn chỉnh phù hợp
  • Mỗi biểu đồ sử dụng theme_cowplot() và một biểu đồ có trục y được di chuyển sang phía bên phải của biểu đồ
#load package
pacman::p_load(cowplot)

# Make first plot of epicurve histogram
#######################################
plot_cases <- linelist %>% 
  
  # plot cases per week
  ggplot()+
  
  # create histogram  
  geom_histogram(
    
    mapping = aes(x = date_onset),
    
    # bin breaks every week beginning monday before first case, going to monday after last case
    breaks = weekly_breaks_all)+  # pre-defined vector of weekly dates (see top of ggplot section)
        
  # specify beginning and end of date axis to align with other plot
  scale_x_date(
    limits = c(min(weekly_breaks_all), max(weekly_breaks_all)))+  # min/max of the pre-defined weekly breaks of histogram
  
  # labels
  labs(
      y = "Daily cases",
      x = "Date of symptom onset"
    )+
  theme_cowplot()


# make second plot of percent died per week
###########################################
plot_deaths <- linelist %>%                        # begin with linelist
  group_by(week = floor_date(date_onset, "week")) %>%  # create week column
  
  # summarise to get weekly percent of cases who died
  summarise(n_cases = n(),
            died = sum(outcome == "Death", na.rm=T),
            pct_died = 100*died/n_cases) %>% 
  
  # begin plot
  ggplot()+
  
  # line of weekly percent who died
  geom_line(                                # create line of percent died
    mapping = aes(x = week, y = pct_died),  # specify y-height as pct_died column
    stat = "identity",                      # set line height to the value in pct_death column, not the number of rows (which is default)
    size = 2,
    color = "black")+
  
  # Same date-axis limits as the other plot - perfect alignment
  scale_x_date(
    limits = c(min(weekly_breaks_all), max(weekly_breaks_all)))+  # min/max of the pre-defined weekly breaks of histogram
  
  
  # y-axis adjustments
  scale_y_continuous(                # adjust y-axis
    breaks = seq(0,100, 10),         # set break intervals of percent axis
    limits = c(0, 100),              # set extent of percent axis
    position = "right")+             # move percent axis to the right
  
  # Y-axis label, no x-axis label
  labs(x = "",
       y = "Percent deceased")+      # percent axis label
  
  theme_cowplot()                   # add this to make the two plots merge together nicely

Bây giờ sử dụng cowplot để chồng lên hai biểu đồ. Sự chú ý đã được tập trung vào việc căn chỉnh trục x, cạnh của trục y và việc sử dụng theme_cowplot().

aligned_plots <- cowplot::align_plots(plot_cases, plot_deaths, align="hv", axis="tblr")
ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])

32.7 Số mới mắc tích lũy

Lưu ý: Nếu sử dụng incidence2, hãy xem chương về cách bạn có thể tính số mới mắc tích lũy bằng một hàm đơn giản. Chương này sẽ đề cập đến cách tính số mới mắc tích lũy và vẽ biểu đồ bằng ggplot().

Nếu bắt đầu bằng bộ số liệu linelist, hãy tạo một cột mới chứa số ca bệnh tích lũy mỗi ngày trong đợt bùng phát bằng cách sử dụng cumsum() từ base R:

cumulative_case_counts <- linelist %>% 
  count(date_onset) %>%                # count of rows per day (returned in column "n")   
  mutate(                         
    cumulative_cases = cumsum(n)       # new column of the cumulative number of rows at each date
    )

10 hàng đầu tiên được hiển thị bên dưới:

Sau đó, cột tích lũy này có thể được vẽ dựa trên date_onset, sử dụng geom_line():

plot_cumulative <- ggplot()+
  geom_line(
    data = cumulative_case_counts,
    aes(x = date_onset, y = cumulative_cases),
    size = 2,
    color = "blue")

plot_cumulative

Nó cũng có thể được phủ đè lên đường cong dịch bệnh, với trục kép bằng cách sử dụng phương pháp cowplot được mô tả ở trên và trong chương Các tips với ggplot:

#load package
pacman::p_load(cowplot)

# Make first plot of epicurve histogram
plot_cases <- ggplot()+
  geom_histogram(          
    data = linelist,
    aes(x = date_onset),
    binwidth = 1)+
  labs(
    y = "Daily cases",
    x = "Date of symptom onset"
  )+
  theme_cowplot()

# make second plot of cumulative cases line
plot_cumulative <- ggplot()+
  geom_line(
    data = cumulative_case_counts,
    aes(x = date_onset, y = cumulative_cases),
    size = 2,
    color = "blue")+
  scale_y_continuous(
    position = "right")+
  labs(x = "",
       y = "Cumulative cases")+
  theme_cowplot()+
  theme(
    axis.line.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    axis.ticks = element_blank())

Bây giờ sử dụng cowplot để chồng lên hai biểu đồ. Sự chú ý đã được tập trung vào việc căn chỉnh trục x, cạnh của trục y và việc sử dụng theme_cowplot().

aligned_plots <- cowplot::align_plots(plot_cases, plot_cumulative, align="hv", axis="tblr")
ggdraw(aligned_plots[[1]]) + draw_plot(aligned_plots[[2]])

32.8 Nguồn tham khảo