Gọi package

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

pacman::p_load(
  rio,          # importing data  
  here,         # relative file pathways  
  janitor,      # data cleaning and tables
  lubridate,    # working with dates
  epikit,       # age_categories() function
  apyramid,     # age pyramids
  tidyverse,    # data manipulation and visualization
  RColorBrewer, # color palettes
  formattable,  # fancy tables
  kableExtra    # table formatting
)

Nhập dữ liệu

Chúng ta sẽ nhập bộ dữ liệu mẫu về những người tiếp xúc và thông tin “theo dõi” của họ. Những dữ liệu này đã được truy xuất và không được lồng ghép với Go.Data API, đồng thời được lưu trữ dưới dạng tệp “.rds”.

Bạn có thể tải xuống tất cả dữ liệu mẫu cho sổ tay này từ chương Tải sách và dữ liệu.

Nếu bạn muốn tải xuống dữ liệu truy vết tiếp xúc mẫu cụ thể cho chương này, hãy sử dụng ba liên kết tải xuống dưới đây:

Bấm để tải xuống dữ liệu điều tra ca mắc (file .rds)

Bấm để tải xuống dữ liệu khai báo tiếp xúc (file .rds)

Bấm để tải xuống dữ liệu theo dõi tiếp xúc (file .rds)

Ở dạng gốc các các tệp có thể tải xuống được, dữ liệu thể hiện dữ liệu gốc được cung cấp bởi Go.Data API (tìm hiểu thêm về APIs ở đây). Đối với mục đích ví dụ ở đây, chúng ta sẽ làm sạch dữ liệu để giúp bạn đọc chương này dễ dàng hơn. Nếu bạn đang sử dụng một phiên bản Go.Data, bạn có thể xem hướng dẫn đầy đủ về cách truy xuất dữ liệu của mình tại đây.

Dưới đây, các bộ dữ liệu được nhập bằng hàm import() từ package rio. Xem chương về Nhập xuất dữ liệu để biết những cách nhập dữ liệu khác nhau. Chúng tôi sử dụng hàm here() để xác định đường dẫn tệp - bạn nên cung cấp đường dẫn tệp cụ thể cho máy tính của mình. Sau đó, chúng tôi sử dụng hàm select() để chỉ chọn một số cột nhất định của dữ liệu, nhằm đơn giản hóa cho mục đích diễn giải.

cases <- import(here("data", "godata", "cases_clean.rds")) %>% 
  select(case_id, firstName, lastName, gender, age, age_class,
         occupation, classification, was_contact, hospitalization_typeid)

contacts <- import(here("data", "godata", "contacts_clean.rds")) %>% 
  mutate(age_class = forcats::fct_rev(age_class)) %>% 
  select(contact_id, contact_status, firstName, lastName, gender, age,
         age_class, occupation, date_of_reporting, date_of_data_entry,
         date_of_last_exposure = date_of_last_contact,
         date_of_followup_start, date_of_followup_end, risk_level, was_case, admin_2_name) %>% 
  mutate(admin_2_name = replace_na(admin_2_name, "Djembe"))

followups <- rio::import(here::here("data", "godata", "followups_clean.rds")) %>% 
  select(contact_id, followup_status, followup_number,
         date_of_followup, admin_2_name, admin_1_name) %>% 
  mutate(followup_status = str_to_lower(followup_status))

relationships <- rio::import(here::here("data", "godata", "relationships_clean.rds")) %>% 
  select(source_visualid, source_gender, source_age, date_of_last_contact,
         date_of_data_entry, target_visualid, target_gender,
         target_age, exposure_type)

Nhân khẩu học

Như được trình bày trong chương Tháp dân số và thang đo Likert, bạn có thể trực quan hóa phân bố theo độ tuổi và giới tính (ở đây chúng tôi sử dụng package apyramid).

Tuổi và Giới tính của người tiếp xúc

Kim tự tháp dưới đây so sánh sự phân bố độ tuổi của những người tiếp xúc theo giới tính. Lưu ý rằng những người tiếp xúc bị missing tuổi được đưa vào thanh riêng của họ ở trên cùng. Bạn có thể thay đổi hiển thị mặc định này, nhưng sau đó hãy xem xét liệt kê số bị missing trong chú thích.

apyramid::age_pyramid(
  data = contacts,                                   # use contacts dataset
  age_group = "age_class",                           # categorical age column
  split_by = "gender") +                             # gender for halfs of pyramid
  labs(
    fill = "Gender",                                 # title of legend
    title = "Age/Sex Pyramid of COVID-19 contacts")+ # title of the plot
  theme_minimal()                                    # simple background

Với cấu trúc dữ liệu Go.Data, bộ dữ liệu relationships bao gồm độ tuổi của cả các ca mắc và những người tiếp xúc, vì vậy bạn có thể sử dụng bộ dữ liệu đó và tạo tháp tuổi cho thấy sự khác biệt giữa hai nhóm người này. Data frame relationships sẽ được biến đổi để chuyển các cột tuổi ở dạng số số thành các nhóm định danh (xem chương Làm sạch số liệu và các hàm quan trọng). Chúng ta cũng sẽ xoay trục dataframe sang định dạng dọc nhằm vẽ biểu đồ dễ dàng hơn với ggplot2 (xem chương Xoay trục dữ liệu).

relation_age <- relationships %>% 
  select(source_age, target_age) %>% 
  transmute(                              # transmute is like mutate() but removes all other columns not mentioned
    source_age_class = epikit::age_categories(source_age, breakers = seq(0, 80, 5)),
    target_age_class = epikit::age_categories(target_age, breakers = seq(0, 80, 5)),
    ) %>% 
  pivot_longer(cols = contains("class"), names_to = "category", values_to = "age_class")  # pivot longer


relation_age

## # A tibble: 200 x 2
##    category         age_class
##    <chr>            <fct>    
##  1 source_age_class 80+      
##  2 target_age_class 15-19    
##  3 source_age_class <NA>     
##  4 target_age_class 50-54    
##  5 source_age_class <NA>     
##  6 target_age_class 20-24    
##  7 source_age_class 30-34    
##  8 target_age_class 45-49    
##  9 source_age_class 40-44    
## 10 target_age_class 30-34    
## # ... with 190 more rows

Bây giờ chúng ta có thể vẽ biểu đồ bộ dữ liệu đã chuyển đổi này với hàm age_pyramid() như trước đây, nhưng thay thế gender bởi category (người tiếp xúc, hoặc ca mắc).

apyramid::age_pyramid(
  data = relation_age,                               # use modified relationship dataset
  age_group = "age_class",                           # categorical age column
  split_by = "category") +                           # by cases and contacts
  scale_fill_manual(
    values = c("orange", "purple"),                  # to specify colors AND labels
    labels = c("Case", "Contact"))+
  labs(
    fill = "Legend",                                           # title of legend
    title = "Age/Sex Pyramid of COVID-19 contacts and cases")+ # title of the plot
  theme_minimal()                                              # simple background

Chúng ta cũng có thể xem xét các đặc điểm khác như cơ cấu nghề nghiệp (ví dụ: ở dạng biểu đồ tròn).

# Clean dataset and get counts by occupation
occ_plot_data <- cases %>% 
  mutate(occupation = forcats::fct_explicit_na(occupation),  # make NA missing values a category
         occupation = forcats::fct_infreq(occupation)) %>%   # order factor levels in order of frequency
  count(occupation)                                          # get counts by occupation
  
# Make pie chart
ggplot(data = occ_plot_data, mapping = aes(x = "", y = n, fill = occupation))+
  geom_bar(width = 1, stat = "identity") +
  coord_polar("y", start = 0) +
  labs(
    fill = "Occupation",
    title = "Known occupations of COVID-19 cases")+
  theme_minimal() +                    
  theme(axis.line = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank())

Số người tiếp xúc trên mỗi ca mắc

Số người tiếp xúc trên mỗi ca mắc có thể là một thước đo quan trọng để đánh giá chất lượng của việc điều tra số người tiếp xúc và mức độ tuân thủ của người dân với việc ứng phó sức khỏe cộng đồng.

Tùy thuộc vào cấu trúc dữ liệu của bạn, điều này có thể được đánh giá thông qua bộ dữ liệu chứa thông tin tất cả các ca mắc và những người tiếp xúc. Trong bộ dữ liệu Go.Data, liên kết giữa các ca mắc (“nguồn”) và những người tiếp xúc (“mục tiêu”) được lưu trữ trong bộ dữ liệu relationships.

Trong bộ dữ liệu này, mỗi hàng là một người tiếp xúc và trong đó có liệt kê ca mắc nguồn. Không có người tiếp xúc nào có mối quan hệ với nhiều ca mắc, nhưng nếu điều này tồn tại, bạn có thể cần tính toán những người đó trước khi vẽ biểu đồ (và khám phá chúng nữa!).

Chúng ta bắt đầu bằng cách đếm số hàng (người tiếp xúc) trên mỗi ca mắc nguồn. Kết quả này được lưu dưới dạng một data frame.

contacts_per_case <- relationships %>% 
  count(source_visualid)

contacts_per_case

## # A tibble: 23 x 2
##    source_visualid     n
##    <chr>           <int>
##  1 CASE-2020-0001     13
##  2 CASE-2020-0002      5
##  3 CASE-2020-0003      2
##  4 CASE-2020-0004      4
##  5 CASE-2020-0005      5
##  6 CASE-2020-0006      3
##  7 CASE-2020-0008      3
##  8 CASE-2020-0009      3
##  9 CASE-2020-0010      3
## 10 CASE-2020-0012      3
## # ... with 13 more rows

Chúng ta sử dụng hàm geom_histogram() để vẽ các dữ liệu này dưới dạng biểu đồ histogram.

ggplot(data = contacts_per_case)+        # begin with count data frame created above
  geom_histogram(mapping = aes(x = n))+  # print histogram of number of contacts per case
  scale_y_continuous(expand = c(0,0))+   # remove excess space below 0 on y-axis
  theme_light()+                         # simplify background
  labs(
    title = "Number of contacts per case",
    y = "Cases",
    x = "Contacts per case"
  )

Làm sạch dữ liệu

Những dữ liệu này có thể tồn tại ở nhiều định dạng khác nhau. Chúng có thể tồn tại dưới dạng trang tính Excel định dạng “ngang” với một hàng cho mỗi người tiếp xúc và một cột cho mỗi “ngày” theo dõi. Xem chương Xoay trục dữ liệu để hiểu về mô tả dữ liệu “dọc” và “ngang”, và cả cách xoay trục dữ liệu sang định dạng ngang hoặc dọc.

Trong ví dụ Go.Data của chúng tôi, những dữ liệu này được lưu trữ trong data frame followups, data frame này có định dạng “dọc” với một hàng cho mỗi tương tác theo dõi. 50 hàng đầu tiên như sau:

CẨN TRỌNG: Cẩn thận với các bản trùng lặp khi xử lý dữ liệu theo dõi; vì có thể có một vài lần theo dõi sai sót trong cùng một ngày cho một người tiếp xúc nhất định. Nó dường như có vẻ là một lỗi nhưng lại phản ánh đúng thực tế - ví dụ: người theo dõi trường hợp tiếp xúc có thể gửi biểu mẫu theo dõi từ sớm trong ngày dù họ chưa liên hệ được với người tiếp xúc và gửi biểu mẫu thứ hai khi họ đã liên hệ được sau đó. Việc này sẽ phụ thuộc vào quy trình thực hiện hoạt động đối với cách bạn muốn xử lý các bản trùng lặp - chỉ cần đảm bảo ghi lại cách tiếp cận của bạn một cách rõ ràng.

Hãy xem chúng ta có bao nhiêu trường hợp hàng “trùng lặp”:

followups %>% 
  count(contact_id, date_of_followup) %>%   # get unique contact_days
  filter(n > 1)                             # view records where count is more than 1  

## # A tibble: 3 x 3
##   contact_id date_of_followup     n
##   <chr>      <date>           <int>
## 1 <NA>       2020-09-03           2
## 2 <NA>       2020-09-04           2
## 3 <NA>       2020-09-05           2

Trong dữ liệu ví dụ của chúng ta, các bản ghi duy nhất áp dụng điều này là những bản ghi missing ID! Chúng ta có thể loại bỏ chúng. Tuy nhiên, với mục đích diễn giải, chúng ta sẽ trình bày các bước nhằm loại bỏ trùng lặp để mỗi người tiếp xúc chỉ có một lần theo dõi truy vết mỗi ngày. Xem chương Loại bỏ trùng lặp để biết thêm chi tiết. Chúng ta sẽ giả định rằng bản ghi cuộc gọi truy vết mới nhất là bản ghi chính xác. Chúng ta cũng tận dụng cơ hội để làm sạch cột followup_number (“ngày” theo dõi sẽ nằm trong khoảng 1 - 14).

followups_clean <- followups %>%
  
  # De-duplicate
  group_by(contact_id, date_of_followup) %>%        # group rows per contact-day
  arrange(contact_id, desc(date_of_followup)) %>%   # arrange rows, per contact-day, by date of follow-up (most recent at top)
  slice_head() %>%                                  # keep only the first row per unique contact id  
  ungroup() %>% 
  
  # Other cleaning
  mutate(followup_number = replace(followup_number, followup_number > 14, NA)) %>% # clean erroneous data
  drop_na(contact_id)                               # remove rows with missing contact_id

Đối với mỗi cuộc gọi truy vết tiếp theo, chúng ta có một trạng thái theo dõi (chẳng hạn như cuộc gọi truy vết có xảy ra hay không và nếu có, người tiếp xúc có triệu chứng hay không). Để xem tất cả các giá trị, chúng ta có thể chạy nhanh lệnh tabyl() (từ janitor) hoặc table() (từ base R) (xem chương Bảng mô tả) bằng followup_status để xem tần suất của từng kết quả.

Trong bộ dữ liệu này, “seen_not_ok” có nghĩa là “được thấy có triệu chứng” và “seen_ok” có nghĩa là “được thấy không có triệu chứng”.

followups_clean %>% 
  tabyl(followup_status)

##  followup_status   n    percent
##           missed  10 0.02325581
##    not_attempted   5 0.01162791
##    not_performed 319 0.74186047
##      seen_not_ok   6 0.01395349
##          seen_ok  90 0.20930233

Vẽ biểu đồ theo thời gian

Vì dữ liệu ngày là liên tục, chúng ta sẽ sử dụng biểu đồ histogram để vẽ chúng với biến date_of_followup được gán cho trục x. Chúng ta có thể vẽ được biểu đồ histogram “xếp chồng” bằng cách chỉ định đối số fill = trong aes(), đối số mà chúng ta gán cho cột followup_status. Do đó, bạn có thể thiết lập chú thích bằng cách sử dụng đối số fill = của labs().

Chúng ta có thể thấy rằng những người tiếp xúc được xác định theo từng đợt (có lẽ là tương ứng với các đợt dịch) và việc hoàn thành theo dõi dường như không cải thiện trong suốt đợt dịch.

ggplot(data = followups_clean)+
  geom_histogram(mapping = aes(x = date_of_followup, fill = followup_status)) +
  scale_fill_discrete(drop = FALSE)+   # show all factor levels (followup_status) in the legend, even those not used
  theme_classic() +
  labs(
    x = "",
    y = "Number of contacts",
    title = "Daily Contact Followup Status",
    fill = "Followup Status",
    subtitle = str_glue("Data as of {max(followups$date_of_followup, na.rm=T)}"))   # dynamic subtitle

CẨN TRỌNG: Nếu bạn đang chuẩn bị nhiều biểu đồ (ví dụ: cho nhiều khu vực pháp lý), bạn sẽ muốn các chú thích xuất hiện giống hệt nhau ngay cả với các mức độ hoàn thiện dữ liệu hoặc thành phần dữ liệu khác nhau. Có thể có những biểu đồ mà không phải tất cả các trạng thái theo dõi đều có trong dữ liệu, nhưng bạn vẫn muốn các danh mục đó xuất hiện trong chú thích. Trong ggplots (như trên), bạn có thể chỉ định đối số drop = FALSE của hàm scale_fill_discrete(). Trong bảng, hãy sử dụng tabyl() để hiển thị số lượng cho tất cả các thực bậc của factor hoặc nếu sử dụng count() từ dplyr, hãy thêm đối số .drop = FALSE để bao gồm số lượng của tất cả các thứ bậc factors.

Theo dõi cá nhân hàng ngày

Nếu sự bùng phát dịch của bạn đủ nhỏ, bạn có thể muốn xem xét từng người tiếp xúc và xem trạng thái của họ trong suốt quá trình theo dõi. May mắn rằng, bộ dữ liệu followups này đã chứa sẵn một cột với “số” ngày theo dõi (1-14). Nếu cột này không tồn tại trong dữ liệu của bạn, bạn có thể tạo ra nó bằng cách tính toán sự khác biệt giữa ngày truy vết và ngày dự định bắt đầu theo dõi người tiếp xúc.

Một cơ chế trực quan hóa thuận tiện (nếu số lượng ca mắc không quá lớn) có thể là một biểu đồ nhiệt, được tạo bằng geom_tile(). Xem thêm chi tiết trong chương Biểu đồ nhiệt.

ggplot(data = followups_clean)+
  geom_tile(mapping = aes(x = followup_number, y = contact_id, fill = followup_status),
            color = "grey")+       # grey gridlines
  scale_fill_manual( values = c("yellow", "grey", "orange", "darkred", "darkgreen"))+
  theme_minimal()+
  scale_x_continuous(breaks = seq(from = 1, to = 14, by = 1))

Phân tích theo nhóm

Có lẽ những dữ liệu theo dõi này đang được kiểm tra hàng ngày hoặc hàng tuần để đưa ra quyết định hành động. Bạn có thể muốn những phân tách có ý nghĩa hơn theo khu vực địa lý hoặc theo nhóm truy vết tiếp xúc. Chúng ta có thể làm điều này bằng cách điều chỉnh các cột được cung cấp tới hàm group_by().

plot_by_region <- followups_clean %>%                                        # begin with follow-up dataset
  count(admin_1_name, admin_2_name, followup_status) %>%   # get counts by unique region-status (creates column 'n' with counts)
  
  # begin ggplot()
  ggplot(                                         # begin ggplot
    mapping = aes(x = reorder(admin_2_name, n),     # reorder admin factor levels by the numeric values in column 'n'
                  y = n,                            # heights of bar from column 'n'
                  fill = followup_status,           # color stacked bars by their status
                  label = n))+                      # to pass to geom_label()              
  geom_col()+                                     # stacked bars, mapping inherited from above 
  geom_text(                                      # add text, mapping inherited from above
    size = 3,                                         
    position = position_stack(vjust = 0.5), 
    color = "white",           
    check_overlap = TRUE,
    fontface = "bold")+
  coord_flip()+
  labs(
    x = "",
    y = "Number of contacts",
    title = "Contact Followup Status, by Region",
    fill = "Followup Status",
    subtitle = str_glue("Data as of {max(followups_clean$date_of_followup, na.rm=T)}")) +
  theme_classic()+                                                                      # Simplify background
  facet_wrap(~admin_1_name, strip.position = "right", scales = "free_y", ncol = 1)      # introduce facets 

plot_by_region