25 接触者の追跡

この章では、実地疫学調査で収集された追跡データの記述統計を行い、接触者追跡データ(以下、追跡データ)を分析する際に考慮すべき重要事項、ならびに課題に直面した時の対処法を説明します。

この章で扱う内容は、他の章で取り扱われているデータ管理と可視化の基礎(データ前処理や整備、表の作り方や、時系列分析など)を参照していますが、その中でも特に、実地疫学調査業務における意思決定に役立つ追跡データに特有の例を取り上げます。例えば、追跡データを時系列や地域ごとに可視化する図や、調査の責任者に報告するための重要業績評価指標(Key Performance Indicator: KPI)の表などを作成します。

本章では、Go.Data プロジェクトのウェブサイトからダウンロードできる追跡データを例として使用しますが、本章で扱う内容は、他のウェブサイトから取得した追跡データを使用しても再現可能です。その場合は、データの構造によって、必要な前処理が異なる可能性があります。

Go.Data の詳細については、Go.Data 公式ドキュメント公式ウェブサイト をご覧ください。

25.1 データ準備

パッケージの読み込み

以下のコードを実行すると、分析に必要なパッケージが読み込まれます。このハンドブックでは、パッケージを読み込むために、pacman パッケージの p_load() を主に使用しています。p_load() は、必要に応じてパッケージをインストールし、現在の R セッションで使用するためにパッケージを読み込む関数です。また、すでにインストールされたパッケージは、R の基本パッケージである base (以下、base R)の library() を使用して読み込むこともできます。R のパッケージに関する詳細は R の基礎 の章をご覧ください。

pacman::p_load(
  rio,          # データのインポート  
  here,         # 相対パスの設定  
  janitor,      # データクリーニング
  lubridate,    # 日付型データ
  epikit,       # age_categories() 
  apyramid,     # 人口ピラミッド
  tidyverse,    # データの処理と可視化
  RColorBrewer, # 配色デザイン
  formattable,  # 見やすい表の作成
  kableExtra    # 表のフォーマットの整備
)

データのインポート

まず、接触者と接触者追跡状況(フォローアップ)のサンプルデータセットをインポートします。この章で使用するサンプルデータは、ネストされていない形で Go.Data API から取得され、“.rds” ファイル形式です。

ハンドブックとデータのダウンロード の章から、このハンドブックで使用するすべてのデータをダウンロードできます。

この章で使用するサンプルデータのみをダウンロードする場合は、以下のリンクをご利用ください。

クリックしてダウンロード 感染症例サンプルデータ(.rds ファイル)

クリックしてダウンロード 接触者サンプルデータ (.rds ファイル)

クリックしてダウンロード 追跡サンプルデータ (.rds ファイル)

上のリンクからダウンロードしたファイルは、Go.Data API(API の詳細についてはこちら)が提供するデータであり、加工されていないため、データをインポートした後、この章で使いやすいようにデータの前処理を行います。Go.Data データの取得について、詳細が知りたい方は、こちら をご覧ください。

以下のコードでは、rio パッケージの import() を使用してサンプルデータセットをインポートしていますが、データをインポートする方法は他にも数多くあります。詳細を知りたい方は、 インポートとエクスポート の章を参照してください。以下のコードでは、here() を使用してファイルパスを指定していますが、コードを実行する際は、使用しているコンピュータ固有のファイルパスに変更してください。データをインポートした後、select() を使用して特定の列のみを選択し、必要のない変数を除外します。

感染症例データ

このサンプルデータ(cases データセット)には、感染者の情報が含まれています。

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

以下は、このデータセットに含まれている nrow(cases) 人の感染者登録データの一覧です。

接触者データ

このサンプルデータ(contacts データセット)は、接触者に関わる情報が含まれた表です。以下のコードを実行する際も、お手元の環境に適したファイルパスを指定してください。データをインポート後、以下の順に従ってデータクリーニングを行います。

  • age_class 変数を因子(ファクタ)型に設定し、レベルを逆にならべ、若い年齢が最初に来るようにする
  • 特定の列のみを選択し、選択した列のうち1つの列の名前を変更する
  • admin_2_name 変数で空欄があった場合、 “Djembe” を割り当て、後ほど作成するデータ可視化の例をわかりやすくする
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"))

以下は、このデータセットに含まれている nrow(contacts) 人の接触者データの一覧です。

追跡データ

このサンプルデータ(followups データセット)には、「フォローアップ」の記録が含まれています。各接触者は、感染への曝露後 14 日間、毎日 1 回、連絡を取る必要があります。

データをインポートした後、いくつかのデータ加工を行います。特定の列を選択し、文字列変数である followup_status をすべて小文字に変換します。

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))

以下では、追跡データの最初の nrow(followups) 人を表示しています。 (各行は 追跡調査の 1 回の記録を表し、追跡調査の結果は 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)

以下は、「関連データ」の最初の 50 行です。感染症例と接触の関係がすべて含まれています。

25.2 記述統計

このハンドブックの他の章で説明されている分析手法や R コードを使用し、感染者、接触者、そして感染者と接触者の関連について記述的な分析を行うことができます。以下に、いくつか例を示します。

人口統計

人口ピラミッドとリッカート尺度 の章で紹介したように、年齢や性別の分布を可視化することができます(ここでは、apyramid パッケージを使用しています)。

接触者の年齢と性別

以下の人口ピラミッドは、接触者の年齢分布を男女別に比較したものです。年齢が不明の接触者は、一番上の unknown に含まれていることに注意してください。年齢不明の接触者を人口ピラミッドから除外することもできますが、その場合は、何人除外されたのかを、プロット下部に注意書きとして記すことをおすすめします。

apyramid::age_pyramid(
  data = contacts,                                   # 接触者サンプルデータを使用
  age_group = "age_class",                           # 年齢変数(因子型)を指定
  split_by = "gender") +                             # 性別による比較
  labs(
    fill = "Gender",                                 # 凡例のタイトル
    title = "Age/Sex Pyramid of COVID-19 contacts")+ # 図のタイトル
  theme_minimal()                                    # 背景テーマの設定

他にも、感染者と接触者両方の年齢が含まれている Go.Data の関連データ(relationships データセット)を使用すると、感染者と接触者の年齢層の違いを表す人口ピラミッドを作成することもできます。関連データを用いて人口ピラミッドを作成する場合は、年齢変数をカテゴリー化し、数字型から因子(ファクタ)型にする必要があります(詳しくは、データクリーニングと主要関数の章を参照ください)。また、ggplot2 パッケージで図をプロットしやすくするために、データを縦型(long型)に変換する必要があります(詳細は、データの縦横変換をご覧ください)。

relation_age <- relationships %>% 
  select(source_age, target_age) %>% 
  transmute(                              # transmute() は mutate() と基本は同じだが、言及されていないすべての列を排除する機能を含む
    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")  # データを縦型(long 型)に変換する


relation_age
## # A tibble: 200 × 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    
## # ℹ 190 more rows

前処理した関連データ(relationships データセット)を使用し、先ほどと同じように age_pyramid() を使用して人口ピラミッドをプロットしてみましょう。ただし、 gender 変数を使用するのではなく、 category 変数(感染者か接触者のどちらであるかを示す変数)を使用する必要があります。

apyramid::age_pyramid(
  data = relation_age,                               # 上のコードで作成したデータセットを使用
  age_group = "age_class",                           # 年齢変数(因子型)を指定
  split_by = "category") +                           # 感染者と接触者でグループ分け
  scale_fill_manual(
    values = c("orange", "purple"),                  # 各グループの色と名前を指定
    labels = c("Case", "Contact"))+
  labs(
    fill = "Legend",                                           # 凡例のタイトル
    title = "Age Pyramid of COVID-19 contacts and cases")+ # 図のタイトル
  theme_minimal()                                              # 背景テーマの設定

他にも、感染者の職業の内訳などの特徴も可視化することができます(ここでは、感染症例データ(cases のデータセット)を使用し、円グラフを作成します)。

# データを前処理し、職業ごとの感染者数を産出する
occ_plot_data <- cases %>% 
  mutate(occupation = forcats::fct_explicit_na(occupation),  # 欠損値を可視化する(NA をカテゴリーとする)
         occupation = forcats::fct_infreq(occupation)) %>%   # 頻度順で因子型のレベルを並べかえる
  count(occupation)                                          # 職業ごとの感染者数を算出する
  
# 円グラフを作成する
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())

感染者一人当たりの接触者数

接触者を見つけるための疫学調査の質や、市民がどの程度公衆衛生施策を遵守しているかを評価するための重要な指標として、感染者一人当たりの接触者数があげられます。

データ構造にもよりますが、すべての感染者と接触者を含むデータセットがあれば、感染者一人当たりの接触者数を算出し、評価を行うことができます。Go.Data のデータでは、関連データ(relationships データセット)に、感染者(sources)と接触者(targets)の関連に関する情報が含まれています。

関連データ(relationships データセット)では、各行が接触者の記録を示し、感染元となった感染者についての情報が記載されています。このデータセットには、複数の感染者と関連している接触者はありませんが、もし複数の感染者と関連している接触者がいる場合には、そのような接触者データを精査し、事前に処理を行う必要があります(この章では触れていません)。

まず、感染者一人当たりの接触者数を算出し、データフレーム(data frame)として保存します。

contacts_per_case <- relationships %>% 
  count(source_visualid)

contacts_per_case
##    source_visualid  n
## 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
## 11  CASE-2020-0013  5
## 12  CASE-2020-0014  3
## 13  CASE-2020-0016  3
## 14  CASE-2020-0018  4
## 15  CASE-2020-0022  3
## 16  CASE-2020-0023  4
## 17  CASE-2020-0030  3
## 18  CASE-2020-0031  3
## 19  CASE-2020-0034  4
## 20  CASE-2020-0036  1
## 21  CASE-2020-0037  3
## 22  CASE-2020-0045  3
## 23            <NA> 17

geom_histogram() を使用し、ヒストグラムを作成します。

ggplot(data = contacts_per_case)+        # 上のコードで作成したデータセットを使用
  geom_histogram(mapping = aes(x = n))+  # 感染者一人当たりの接触者数のヒストグラムを作成
  scale_y_continuous(expand = c(0,0))+   # y軸の0より下の余分なスペースを削除
  theme_light()+                         # 図の背景テーマの設定
  labs(
    title = "Number of contacts per case",
    y = "Cases",
    x = "Contacts per case"
  )

25.3 接触者追跡調査(フォローアップ)

ほとんどの追跡データには、隔離された人の毎日の症状チェックの結果を記録した接触者調査(フォローアップ)のデータが含まれています。接触者調査の記録を分析することで、今後の公衆衛生施策の策定に役立ち、追跡できなくなる危険性のある接触者や、疾患を発症する危険性のある接触者を特定することができます。

データの前処理

接触者調査に関する情報は、様々なフォーマットで記録されています。例えば、Excel シートに各接触者を 1 行ごとに記録し、調査の結果を1日ごとに列に記録した、横型(wide 型)データなどです。縦型(long 型)データと横型(wide 型)の詳細や、データの縦横変換の方法については、 データの縦横変換 の章を参照してください。

Go.Data からダウンロードしたデータは、追跡データ(followups データセット)に含まれており、調査の記録が各行に記載されている縦型データとなっています。以下では、接触者追跡調査データの最初の50行を表示しています。

<style=“color:orange;”>注意: 接触者調査(フォローアップ)に関するデータを扱う際は、重複した調査記録に注意してください。ある接触者に対し、同じ日に複数の調査(フォローアップ)が誤って行われる可能性があります。例えば、調査員が接触者と連絡がとれなかった際にその記録を午前中に提出し、後に連絡できた際に 2 つ目の記録を提出することが起こりえます。重複した記録をどのように処理するかについては、調査の運用状況によりますが、重複記録を含んだデータを提示する際は、重複記録をどのように処理したかを明確に記載してください。

追跡データ(followups データセット)に、重複している記録がいくつあるのか、チェックしてみましょう。

followups %>% 
  count(contact_id, date_of_followup) %>%   # 接触者ID毎に、調査日の頻度を算出
  filter(n > 1)                             # 1日2回以上調査が行われた日を表示
##   contact_id date_of_followup n
## 1       <NA>       2020-09-03 2
## 2       <NA>       2020-09-04 2
## 3       <NA>       2020-09-05 2

ここで扱っているサンプルデータでは、ID のない記録だけが該当しています。下のコードで、重複した記録を排除していきます。重複記録を排除することで、1 日に 1 人につき 1 回の調査記録のみ残るようになります。 詳細は、重複データの排除の章をご覧ください。ここでは、最新の調査記録が正しいものであると仮定します。また、下のコードでは、一緒に followup_number 列(調査の「日」ごとに作成され、1 ~ 14 日まである)をクリーニングします。

followups_clean <- followups %>%
  
  # 重複記録の排除
  group_by(contact_id, date_of_followup) %>%        #接触者ID、調査日ごとに記録をグループ化する
  arrange(contact_id, desc(date_of_followup)) %>%   # 接触者ID毎に、調査日の降順で記録を並び替える(日付の新しい順)
  slice_head() %>%                                  # 接触者ID毎に、最新の調査日の記録のみ残す  
  ungroup() %>% 
  
  # 他の前処理
  mutate(followup_number = replace(followup_number, followup_number > 14, NA)) %>% # 調査日のエラーの処理
  drop_na(contact_id)                               # 接触者IDが欠損してる記録を排除する

それぞれの調査記録には、調査の結果(接触者と連絡が取れたか、また、連絡が取れた場合は接触者に症状があったか、など)が記録されています。janitor パッケージの tabyl() 、または base R の table()followup_status 変数を指定して実行すると、調査結果の分布が確認できます(詳細は、記述統計表の作り方の章をご覧ください)。

このデータセットでは、“seen_not_ok” は「いくつか症状がある」、“seen_ok” は「症状がない」という意味です。

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

追跡調査状況を時系列でプロットする

日付データは連続しているため、date_of_followup を x 軸に指定したヒストグラムを使用し、追跡データ(followups のデータセット)を時系列でプロットできます。aes() 内の fill = 引数に followup_status 変数を指定することで、「積み上げ型」のヒストグラムを作成し、labs() 内の fill = 引数を使用して凡例のタイトルを作成することができます。

このデータセットを使用して作成したヒストグラムでは、接触者の分布が波状に確認され(おそらく感染者の流行の波と連動している)、追跡調査が行われていない場合が多く、流行の過程で調査の達成状況が改善されていないことがわかります。

ggplot(data = followups_clean)+
  geom_histogram(mapping = aes(x = date_of_followup, fill = followup_status)) +
  scale_fill_discrete(drop = FALSE)+   # followup_status 変数のすべての値を凡例に表示する(図に表示されていないものも含めて)
  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)}"))   # 動的なサブタイトルをつける
<style=“color: orange;”>注意: 同一コードで複数のプロットを作成する場合(複数の管轄区域に同じプロットを提示する場合など)、データの完成度やデータの構成が異なっていても、凡例は同じように表示されるようにする必要があります。例えば、作成した複数の図の中に followup_status 変数(調査結果変数)のすべてのカテゴリーがデータに含まれていない図があるかもしれませんが、含まれていない結果のカテゴリーについても凡例に表示したい場合です。上のような ggplot を使用した図では、scale_fill_discrete() 内で drop = FALSE に指定すると、図に表示されていないカテゴリーを含むすべてのカテゴリーを凡例に表示することができます。表の場合は、すべてのカテゴリーのカウントを表示する tabyl() を使用するか、または、dplyr パッケージの count().drop = FALSE を追加し、すべてのカテゴリーのカウントを含めることができます。

追跡調査状況を接触者別にプロットする

調査対象のアウトブレイクの規模が小さい場合は、各接触者のフォローアップ状況を個別に確認したい場合があります。この追跡データ(followups データセット)には、調査の日にち「番号」を示す列がすでに含まれています(1 ~ 14 日まで一日ごとに列が作成してあります)。調査日ごとの列がデータにない場合は、曝露日(調査対象の接触者が感染者と接触した日)とその接触者に対して調査を開始した日の差を計算し、1 ~ 14 日まで一日ごとに列を作成する必要があります。

データのわかりやすい視覚化の例として(アウトブレイクの規模が大きくなければ)、geom_tile() を用いたヒートマップがあります。詳細は、ヒートマップ の章をご覧ください。

ggplot(data = followups_clean)+
  geom_tile(mapping = aes(x = followup_number, y = contact_id, fill = followup_status),
            color = "grey")+       # 灰色のグリッド線をひく
  scale_fill_manual( values = c("yellow", "grey", "orange", "darkred", "darkgreen"))+
  theme_minimal()+
  scale_x_continuous(breaks = seq(from = 1, to = 14, by = 1))

追跡調査状況をクループ別にプロットする

このような接触者追跡調査に関するデータは、疫学調査の業務的意思決定のために日ごとまたは週ごとに閲覧されていることが多いと思います。group_by() で指定する列を調整することにより、地域別や調査チーム別など、より意味のある集計結果をプロットすることができます。

plot_by_region <- followups_clean %>%                                        # 前処理した followups のデータセットを使用する
  count(admin_1_name, admin_2_name, followup_status) %>%   # 地域別、調査チーム別のカウントを算出する (新たに 'n' 列を作成する)
  
  # ggplot() を使用して図を作成する
  ggplot(                                         # ggplot
    mapping = aes(x = reorder(admin_2_name, n),     # admin_2_name の因子(ファクタ)レベルを 'n' 列のカウントをもとに並び替える
                  y = n,                            # 'n' 列のカウントをy軸に指定する
                  fill = followup_status,           # フォローアップの結果カテゴリーごとに色付けする
                  label = n))+                      # geom_label() 用              
  geom_col()+                                     # 積み上げ型棒グラフ
  geom_text(                                      # テキストを追加する
    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()+                                                                      # 図の背景のテーマの設定
  facet_wrap(~admin_1_name, strip.position = "right", scales = "free_y", ncol = 1)      # facet_wrap で図の右側に admin_1_name を表示する 

plot_by_region

25.4 重要業績評価指標(KPI)

接触者追跡調査の成果を評価するために、進捗や目標達成状況を様々なレベルで細分化された期間で算出し、モニターする重要業績評価指標(Key Performance Indicator: KPI)があります。KPI には多くの種類がありますが、計算方法と基本的な表のフォーマットを理解すると、様々な KPI を入れたり外したりすることが簡単にできます。

接触者調査の KPI については、ResolveToSaveLives.org など、参考となる資料がたくさんあります。KPI 作成作業の多くは、データ構造を確認し、すべての選択・除外基準を考えることです。以下に、Go.Data を使用した KPI の例を示します。

カテゴリー KPI 分子 分母
プロセスの評価 - 接触者追跡調査の迅速さ 全感染者のうち、感染報告から24時間以内に疫学調査が行われ、隔離された感染者の割合 COUNT OF case_id WHERE (date_of_reporting - date_of_data_entry) < 1 day AND (isolation_startdate - date_of_data_entry) < 1 day COUNT OF case_id
プロセスの評価 - 接触者追跡調査の迅速さ 全接触者のうち、感染者の接触報告から24時間以内に通知され、隔離された接触者の割合 COUNT OF contact_id WHERE followup_status == “SEEN_NOT_OK” OR “SEEN_OK” AND date_of_followup - date_of_reporting < 1 day COUNT OF contact_id
プロセスの評価 - 感染症検査の完了率 全感染者のうち、症状発現から3日以内に検査され、疫学調査が行われた者の割合 COUNT OF case_id WHERE (date_of_reporting - date_of_onset) < =3 days COUNT OF case_id
アウトカムの評価 - 全体評価 全感染者のうち、感染報告前に接触者として報告されていた者の割合 COUNT OF case_id WHERE was_contact == “TRUE” COUNT OF case_id

このセクションでは、接触者追跡調査状況を複数の管轄区域にわたって表示する表を作成していきます。表を作成した後、formattable パッケージを使用し、更に見やすく表を編集していきます(flextable などの他のパッケージを使用することもできます。詳細は、見やすい表の作り方 の章を参照してください)。

表の作成手順は、接触者追跡調査に関するデータの構造によりって変わります。 dplyr パッケージを使用してデータを要約する方法についての詳細は、記述統計表の作り方 の章をご覧ください。

ここでは、データの変化に合わせて連動する表を作成します。結果を面白くするために、特定の日に report_date (この表を提示し、追跡調査の進捗を報告する日)を設定し、表を作成します(この例では、2020年6月10日を選びます)。データは report_date でフィルタリングされ、report_date の日、または、report_date より以前に報告された記録のみ残ります。

# report_date (報告日)を設定する
report_date <- as.Date("2020-06-10")

# report_date (報告日)でデータセットをフィルタリングする
table_data <- followups_clean %>% 
  filter(date_of_followup <= report_date)

report_date でフィルタリング後、以下の手順でデータを整理し、表を作成していきます。

  1. 追跡データ(followups データセット)を選択し、各接触者について、以下 3 つの関心のある指標を算出します。
  • 一番最近行った調査の日付(接触者と連絡が取れたかなどの調査結果は問わない)
  • 接触者と連絡が取れた調査のうち、最も新しい日付
  • 接触者と連絡が取れた最新の調査での、接触者の状況(例:症状あり又は症状なし)
  1. ステップ 1 で作成したデータを、曝露日(感染者と接触した日)などの接触者に関する他の情報を含む接触データに結合します。次に、各接触者について、曝露日からの日数など、関心のある指標を算出します。
  2. 結合したデータを、管轄区域別(admin_2_name)にグループ化し、管轄区域別の要約統計量を算出します。
  3. 最後に、作成した表の形式を整えます。

では、以上の手順を実際に R で実行していきましょう。ステップ 1:追跡データ(followups データセット)を整備し、必要な指標を算出する。

followup_info <- table_data %>% 
  group_by(contact_id) %>% 
  summarise(
    date_last_record   = max(date_of_followup, na.rm=T),
    date_last_seen     = max(date_of_followup[followup_status %in% c("seen_ok", "seen_not_ok")], na.rm=T),
    status_last_record = followup_status[which(date_of_followup == date_last_record)]) %>% 
  ungroup()

作成したデータセットを以下に表示します。

ステップ 2:次に、作成したデータセットを接触者データ(contacts データセット)に結合し、他に必要な指標を算出します。

contacts_info <- followup_info %>% 
  right_join(contacts, by = "contact_id") %>% 
  mutate(
    database_date       = max(date_last_record, na.rm=T),
    days_since_seen     = database_date - date_last_seen,
    days_since_exposure = database_date - date_of_last_exposure
    )

2 つのデータセットを結合すると、データは以下のようになります。接触データに含まれていた列は右側に、算出した関心のある指標は一番右に表示されます。

ステップ 3:結合したデータを、管轄区域別にグループ化し、区域別に要約統計量を算出します。

contacts_table <- contacts_info %>% 
  
  group_by(`Admin 2` = admin_2_name) %>%
  
  summarise(
    `Registered contacts` = n(),
    `Active contacts`     = sum(contact_status == "UNDER_FOLLOW_UP", na.rm=T),
    `In first week`       = sum(days_since_exposure < 8, na.rm=T),
    `In second week`      = sum(days_since_exposure >= 8 & days_since_exposure < 15, na.rm=T),
    `Became case`         = sum(contact_status == "BECAME_CASE", na.rm=T),
    `Lost to follow up`   = sum(days_since_seen >= 3, na.rm=T),
    `Never seen`          = sum(is.na(date_last_seen)),
    `Followed up - signs` = sum(status_last_record == "Seen_not_ok" & date_last_record == database_date, na.rm=T),
    `Followed up - no signs` = sum(status_last_record == "Seen_ok" & date_last_record == database_date, na.rm=T),
    `Not Followed up`     = sum(
      (status_last_record == "NOT_ATTEMPTED" | status_last_record == "NOT_PERFORMED") &
        date_last_record == database_date, na.rm=T)) %>% 
    
  arrange(desc(`Registered contacts`))

ステップ 4:formattableknitr パッケージを使用し、見やすいように表を整えていきます。また、report_date (報告日)に関する脚注を追記します。

contacts_table %>%
  mutate(
    `Admin 2` = formatter("span", style = ~ formattable::style(
      color = ifelse(`Admin 2` == NA, "red", "grey"),
      font.weight = "bold",font.style = "italic"))(`Admin 2`),
    `Followed up - signs`= color_tile("white", "orange")(`Followed up - signs`),
    `Followed up - no signs`= color_tile("white", "#A0E2BD")(`Followed up - no signs`),
    `Became case`= color_tile("white", "grey")(`Became case`),
    `Lost to follow up`= color_tile("white", "grey")(`Lost to follow up`), 
    `Never seen`= color_tile("white", "red")(`Never seen`),
    `Active contacts` = color_tile("white", "#81A4CE")(`Active contacts`)
  ) %>%
  kable("html", escape = F, align =c("l","c","c","c","c","c","c","c","c","c","c")) %>%
  kable_styling("hover", full_width = FALSE) %>%
  add_header_above(c(" " = 3, 
                     "Of contacts currently under follow up" = 5,
                     "Status of last visit" = 3)) %>% 
  kableExtra::footnote(general = str_glue("Data are current to {format(report_date, '%b %d %Y')}"))
Of contacts currently under follow up
Status of last visit
Admin 2 Registered contacts Active contacts In first week In second week Became case Lost to follow up Never seen Followed up - signs Followed up - no signs Not Followed up
Djembe 59 30 44 0 2 15 22 0 0 0
Trumpet 3 1 3 0 0 0 0 0 0 0
Venu 2 0 0 0 2 0 2 0 0 0
Congas 1 0 0 0 1 0 1 0 0 0
Cornet 1 0 1 0 1 0 1 0 0 0
Note:
Data are current to Jun 10 2020

25.5 感染連鎖の可視化

ヒートマップ の章で解説したように、geom_tile() を使用して「誰が誰から感染したか」を可視化するヒートマップを作成することができます。

Go.Data では、接触者データ(contacts データセット)に新しい接触者が追加されると、API によって、関連データ(relationships データセット)に、その接触者と接触した感染者の関係情報が追加されます。relationships データセットには、各接触者と接触した感染者の情報が含まれているので、このデータセットを使用すると比較的簡単にヒートマップを作成することができます。以下に、関連データ(relationships データセット)の最初の 50 行を表示します。

接触者と感染者の年齢を比較した人口ピラミッドを 25.2 記述統計のセクションで作成したように、必要な変数以外を除外し、接触者と感染者の両方について、年齢をカテゴリー化する必要があります。

heatmap_ages <- relationships %>% 
  select(source_age, target_age) %>% 
  mutate(                            # transmuteはmutateと似ていますが、他の列をすべて除外します
    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))) 

まず、ヒートマップ の章で解説したように、クロス集計表を作成します。

cross_tab <- table(
  source_cases = heatmap_ages$source_age_class,
  target_cases = heatmap_ages$target_age_class)

cross_tab
##             target_cases
## source_cases 0-4 5-9 10-14 15-19 20-24 25-29 30-34 35-39 40-44 45-49 50-54
##        0-4     0   0     0     0     0     0     0     0     0     1     0
##        5-9     0   0     1     0     0     0     0     1     0     0     0
##        10-14   0   0     0     0     0     0     0     0     0     0     0
##        15-19   0   0     0     0     0     0     0     0     0     0     0
##        20-24   1   1     0     1     2     0     2     1     0     0     0
##        25-29   1   2     0     0     0     0     0     0     0     0     0
##        30-34   0   0     0     0     0     0     0     0     1     1     0
##        35-39   0   2     0     0     0     0     0     0     0     1     0
##        40-44   0   0     0     0     1     0     2     1     0     3     1
##        45-49   1   2     2     0     0     0     3     0     1     0     3
##        50-54   1   2     1     2     0     0     1     0     0     3     4
##        55-59   0   1     0     0     1     1     2     0     0     0     0
##        60-64   0   0     0     0     0     0     0     0     0     0     0
##        65-69   0   0     0     0     0     0     0     0     0     0     0
##        70-74   0   0     0     0     0     0     0     0     0     0     0
##        75-79   0   0     0     0     0     0     0     0     0     0     0
##        80+     1   0     0     2     1     0     0     0     1     0     0
##             target_cases
## source_cases 55-59 60-64 65-69 70-74 75-79 80+
##        0-4       1     0     0     0     0   0
##        5-9       1     0     0     0     0   0
##        10-14     0     0     0     0     0   0
##        15-19     0     0     0     0     0   0
##        20-24     1     0     0     0     0   1
##        25-29     0     0     0     0     0   0
##        30-34     1     0     0     0     0   0
##        35-39     0     0     0     0     0   0
##        40-44     1     0     0     0     1   1
##        45-49     2     1     0     0     0   1
##        50-54     1     0     1     0     0   1
##        55-59     0     0     0     0     0   0
##        60-64     0     0     0     0     0   0
##        65-69     0     0     0     0     0   0
##        70-74     0     0     0     0     0   0
##        75-79     0     0     0     0     0   0
##        80+       0     0     0     0     0   0

次に、クロス集計表を縦型(long 型)データに変換します。

long_prop <- data.frame(prop.table(cross_tab))

最後に、接触者がどの年齢層の感染者から感染したかを年齢層別に表すヒートマップを作成します。

ggplot(data = long_prop)+       # long型データを使用(Freq変数を割合として使用)
  geom_tile(                    # タイルで可視化
    aes(
      x = target_cases,         # x軸を接触者の年齢層に指定
      y = source_cases,     # y軸を感染者の年齢層に指定
      fill = Freq))+            # Freq変数でタイルを色付け
  scale_fill_gradient(          # タイルの色の調整
    low = "blue",
    high = "orange")+
  theme(axis.text.x = element_text(angle = 90))+
  labs(                         # プロット、軸、凡例のタイトルを指定
    x = "Target case age",
    y = "Source case age",
    title = "Who infected whom",
    subtitle = "Frequency matrix of transmission events",
    fill = "Proportion of all\ntranmsission events"     # 凡例のタイトル
  )