26 標本調査データ分析

26.1 概要

この章では、標本調査データ分析のためのいくつかのパッケージの使用方法を説明します。

ほとんどの調査用 R パッケージは、重み付け分析を行うために surveyパッケージに依存しています。本章では、survey パッケージと同様に srvyr パッケージ（tidyverse スタイルのコーディングを可能にする survey パッケージのラッパー）および gtsummary パッケージ（出版原稿レベルの整った表の出力を可能にする survey パッケージのラッパー）を使用します。オリジナルの survey パッケージは、tidyverse スタイルのコーディングはできませんが、調査重み付け一般化線形モデル（survey-weighted generalised linearmodel）を使用できるという利点があります（後日この章に追加される予定です）。また、sitrep パッケージの関数を使用して、標本抽出法による重み付けを作成するデモも行います（注：このパッケージは現在 CRAN にはありませんが、Github からインストールできます）。

この章のほとんどは、“R4Epis” プロジェクトで行われた作業に基づいています。詳細なコードと R markdown テンプレートファイルについては、“R4Epis” の Github ページを参照してください。survey パッケージを利用したのコードの一部は、EPIET ケーススタディの初期バージョンに基づいています。

現在、この章では標本サイズの計算や標本抽出法そのものについては触れていません。簡単に利用できる標本サイズ計算機については、OpenEpi を参照してください。ハンドブックの GIS の基礎の章には、いずれ空間的無作為抽出のセクションが設けられ、本章には、標本サイズの計算だけでなく、無作為抽出するための「標本抽出枠」（sampling frame）のセクションが設けられる予定です。

調査データ
観察期間
重み付け
調査デザイン対象
記述的分析
重み付けされた割合
重み付けされた比率

26.2 準備

パッケージ

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

## CRAN からパッケージを読み込む
pacman::p_load(rio,          # ファイルのインポート
               here,         # ファイルパスの指定
               tidyverse,    # データ管理 + ggplot2 のグラフ
               tsibble,      # 時系列データセットを扱う
               survey,       # 調査関数
               srvyr,        # survey パッケージ用の dplyr ラッパー
               gtsummary,    # 表作成のための survey パッケージ
               apyramid,     # 年齢ピラミッド作成のためのパッケージ
               patchwork,    # ggplot を組み合わせるためのパッケージ
               ggforce       # 沖積図・サンキーダイアグラム
               ) 

## load packages from github
pacman::p_load_gh(
     "R4EPI/sitrep"          # 観察期間・重み付け関数
)

データの読み込み

本章で使用する例示用データセットです。

架空の死亡率調査データ
調査対象地域の架空の母集団
架空の死亡率調査データのデータ辞書

これらは、MSF OCA 倫理審査委員会が事前に承認した調査に基づいています。架空のデータセットは、“R4Epis” プロジェクトの一環として作成されました。これはすべて、OpenDataKit をベースにしたデータ収集ソフトウェアである KoboToolbox を使って収集したデータに基づいています。

Kobo では、収集したデータだけでなく、そのデータセットのデータ辞書もエクスポートできます。データの前処理が容易になり、変数や質問文を調べるのに便利なので、この方法を強くお勧めします。

ヒント： Kobo のデータ辞書形式（例、MSF-survey-dict.xlsx）では、調査ごとシートに分かれており、各シートの “name” 列に変数名（質問）が記載されています。各変数がとり得る値（回答の選択肢）は、“options” シートで設定されています。“options” シートでは、“name” 列には選択肢の省略形が、“label::english” および “label::french” 列にはそれぞれの言語で適切な長さの選択肢が表示されています。Epidict パッケージの msf_dict_survey() 関数を使って、Kobo データ辞書形式の Excel ファイルをインポートすると、対象の調査シートがデータフレームに変換されるので、調査データをコード記述に簡単に利用できます。

注意：ハンドブックとデータのダウンロードからダウンロードできる例示用データセット（survey_dict.xlsx）は、Kobo からの単純なエクスポートとは出力形式が異なります（Koboでは、異なる質問表をそれぞれ個別にエクスポートします）異なる質問表をマージするには、以下の「調査データ」のセクションを参照してください。

データセットは rio パッケージの import() 関数を使ってインポートされます。データをインポートするさまざまな方法については、データのインポート・エクスポートの章を参照してください。

# 調査データをインポートする
survey_data <- rio::import("survey_data.xlsx")

# データ辞書を R にインポートする
survey_dict <- rio::import("survey_dict.xlsx")

調査データの最初の 10 行を以下に示します。

適切な重みを作成するために、標本を抽出した集団（母集団）についてのデータをインポートしましょう。このデータの形式はいろいろなファイル形式で作成できますが、以下に例示するような形式をお勧めします（このデータは、Excel に直接入力しても作成できます）。

# 母集団データをインポートする
population <- rio::import("population.xlsx")

母集団データの最初の 10 行を以下に示します。

クラスター調査の場合は、クラスター単位で調査の重みを加えることが必要な場合もあります。クラスター単位のデータは上記と同じ手段でインポートできます。あるいは、クラスターや世帯のカウント数が少ない場合は、以下のように入力してクラスターとクラスターが含む世帯数の 2 つの列を定義できます。いずれの方法にしても、調査データと一致するクラスターを識別する列と、各クラスターが含む世帯数を表す列が必要になります。

## 各クラスターの世帯数を定義する
cluster_counts <- tibble(cluster = c("village_1", "village_2", "village_3", "village_4", 
                                     "village_5", "village_6", "village_7", "village_8",
                                     "village_9", "village_10"), 
                         households = c(700, 400, 600, 500, 300, 
                                        800, 700, 400, 500, 500))

データの前処理（データクリーニング）

日付を値として持つ変数（列）が適切な形式であることを以下で確認しています。他にもいくつかの方法がありますが（詳しくは「日付型データ」の章をご覧ください）、データ辞書を使って日付を定義するのが手っ取り早くて簡単です。

また、epikit パッケージの age_categories() 関数を使って、年齢グループの変数を作成します（詳細はデータクリーニングと主要関数の章を参照）。さらに、各クラスターがどの担当保健所に含まれるかを定義する文字型変数（列）を作成します。

最後に、yes/no で回答されている変数（列）をすべて TRUE/FALSE の値を取る変数に上書きします。値が yes/no のままであると、これらの変数は survey パッケージの人口を扱う関数で使用できません。

## survery_dict オブジェクトからtype 列に値 date を持つ行を選択する 
DATEVARS <- survey_dict %>% 
  filter(type == "date") %>% 
  filter(name %in% names(survey_data)) %>% 
  ## 前行でフィルタした行の "name" 列の値が survey_data オブジェクト内の変数名（列名）と、一致した行をフィルタする
  pull(name) # type 列に値 date を持つ行を選択する
  
## survey_data オブジェクト内の変数（列）で DATEVARS に含まれる変数（列）を日付型に変更する 
survey_data <- survey_data %>%
  mutate(across(all_of(DATEVARS), as.Date))


## 年齢が月単位の値（1 歳未満）を年齢の変数に加える（12 で割る）
survey_data <- survey_data %>% 
  mutate(age_years = if_else(is.na(age_years), 
                             age_months / 12, 
                             age_years))

## 年齢層を格納する変数を定義する
survey_data <- survey_data %>% 
     mutate(age_group = age_categories(age_years, 
                                    breakers = c(0, 3, 15, 30, 45)
                                    ))


## 他に cluster_number 変数（列）の内容を基にした health_distinct （担当保健所）という文字型変数（列）を作成する
survey_data <- survey_data %>% 
  mutate(health_district = case_when(
    cluster_number %in% c(1:5) ~ "district_a", 
    TRUE ~ "district_b"
  ))


## survey_dict オブジェクトから type 列 に yn を値として持つ行を選択する 
YNVARS <- survey_dict %>% 
  filter(type == "yn") %>% 
  filter(name %in% names(survey_data)) %>% 
  ## 前行でフィルタした行の "name" 列の値が survey_data オブジェクト内の変数名（列名）と、一致した行をフィルタする
  pull(name) # type 列に値 yn を持つ行を選択する
  
## survey_data オブジェクト内の変数（列）で YNVARS に含まれる変数（列）をロジカル型に変更する 
survey_data <- survey_data %>%
  mutate(across(all_of(YNVARS), 
                str_detect, 
                pattern = "yes"))

26.3 調査データ

標本調査の方法には、さまざまな標本抽出法があります。このセクションでは、以下の標本抽出法のコード実装例を紹介します。 - 層化抽出法 - クラスター抽出法 - 層化クラスター抽出法

「データの読み込み」セクションで例示したように、（調査票をどのように計画するかにもよりますが）各レベル（水準）のデータは、Kobo から別々のデータセットとしてエクスポートされます。下記のコード例では、世帯と、その世帯内の個人がそれぞれ 1 つのレベルで表現されています。

この 2 つのレベルは、一意の識別子で紐付けられています。Kobo のデータセットの場合、識別子は世帯レベルの “_index” という変数の値であり、各個人レベルの “_parent_index” という変数名の値と一致します。left_join() による結合操作で、識別子が一致する個人ごとに世帯の新しい行が作成されます。詳細はハンドブックのデータの結合の章を参照してください。

## 個人データと世帯データを結合し、完全なデータセットを作成する
survey_data <- left_join(survey_data_hh, 
                         survey_data_indiv,
                         by = c("_index" = "_parent_index"))


## 個人と世帯の 2 つのレベルの組み合わせで、新たに識別子を作成する
survey_data <- survey_data %>% 
     mutate(uid = str_glue("{index}_{index_y}"))

26.4 観察期間

死亡率調査において、対象期間における適切な死亡率を算出するために、各個人が対象地域にどのくらいの期間、居住したかという情報が必要です。この居住に関する情報は、すべての調査に当てはまるわけではありませんが、特に死亡率調査の場合は重要です。なぜなら、移動人口や避難民人口を含んだ状況で頻繁に行われる調査だからです。

死亡率調査を行うためには、まず、対象期間を設定する必要があります。対象期間の設定は、想起期間（リコール期間、質問に答える際に参加者にその期間の健康状態を考慮して回答を求める期間）としても知られています。例えば、この対象期間を利用して、対象期間外の死亡報告などは、不適切な日付として欠損値を設定できます。

## 想起期間の開始日・終了日を設定する
## データセットに日付型変数（列）として設定
## (例：転入日や質問日)
survey_data <- survey_data %>% 
  mutate(recall_start = as.Date("2018-01-01"), 
         recall_end   = as.Date("2018-05-01")
  )


# ルールに基づいて不適切な日付に NA を設定する 
## 例：想起期間開始前の転入、想起期間終了後の転出
survey_data <- survey_data %>%
      mutate(
           arrived_date = if_else(arrived_date < recall_start, 
                                 as.Date(NA),
                                  arrived_date),
           birthday_date = if_else(birthday_date < recall_start,
                                  as.Date(NA),
                                  birthday_date),
           left_date = if_else(left_date > recall_end,
                              as.Date(NA),
                               left_date),
           death_date = if_else(death_date > recall_end,
                               as.Date(NA),
                               death_date)
           )

次に、日付型の値を格納している変数（列）を使って、観察期間開始日と観察期間終了日を設定します。sitrep パッケージの find_start_date() 関数と find_end_date() 関数を使って、観察期間開始・終了イベント発生理由とそれぞれの日付を調べます。そして、調べた観察期間開始・終了の日付を使って観察日数（person-time）を計算します。

観察期間開始日：想起期間内の適切なイベント開始の中で最も古い日付。想起期間の開始日（事前に設定する）、または該当するイベントがある場合は想起期間開始以後の日付（例：転入や出生など）が割り当てられる

観察期間終了日：想起期間内の適切なイベント終了の中で最も新しい日付。想起期間の終了日（事前に設定する）、または該当するイベントがある場合は想起期間終了以前の日付（例：転出や死亡など）が割り当てられる

## 観察期間開始日、観察期間終了日、それぞれのイベント発生理由を格納するための新しい変数（列）を作成する
survey_data <- survey_data %>% 
     ## イベント発生理由が、出生、転入（世帯の転入、難民キャンプの転入）のうち
     ## 入力された最も古い日付を選択する
     find_start_date("birthday_date",
                  "arrived_date",
                  period_start = "recall_start",
                  period_end   = "recall_end",
                  datecol      = "startdate",
                  datereason   = "startcause" 
                 ) %>%
     ## イベント発生理由が、転出（難民キャンプの転出ふくむ）、死亡、調査の終了のうち
     ## 入力された最も新しい日付を選択する
     find_end_date("left_date",
                "death_date",
                period_start = "recall_start",
                period_end   = "recall_end",
                datecol      = "enddate",
                datereason   = "endcause" 
               )


## 想起期間開始・終了時点で継続して居住していた記録にラベルをつける（出生・死亡を除く）
survey_data <- survey_data %>% 
     mutate(
       ## 想起期間の開始日を個人の観察期間開始日に設定する（個人の観察期間開始日が欠損している場合への対応） 
       startdate = if_else(is.na(startdate), recall_start, startdate), 
       ## 個人の観察期間開始日が、出生日と等しくなく、かつ、想起期間開始日と等しい場合、
       ## 観察期間開始理由を "Present at start" と設定する 
       startcause = if_else(startdate == recall_start & startcause != "birthday_date",
                              "Present at start", startcause), 
       ## 想起期間の終了日を個人の観察期間終了日に設定する（個人の観察期間終了日が欠損している場合への対応） 
       enddate = if_else(is.na(enddate), recall_end, enddate), 
       ## 個人の観察期間終了日が、死亡日と等しくなく、かつ、想起期間終了日に等しい場合は、
       ## 観察期間終了理由を "Present at end" と設定する
       endcause = if_else(enddate == recall_end & endcause != "death_date", 
                            "Present at end", endcause))


## 観察期間を日単位で計算し obstime 変数（列）に設定する
survey_data <- survey_data %>% 
  mutate(obstime = as.numeric(enddate - startdate))

26.5 調査の重み付け

調査の重み付けを実行する前に、誤った観測値の除外が必要です。例えば、負の観察期間が入力されている観察値がある場合、その詳細を確認する必要があります（sitrep パッケージの assert_positive_timespan() で負の観察期間の記録を一覧できます）。他に必要な前処理は、空行の除外や（例えば drop_na(uid) を使用して）、重複した行の除外です（詳細はハンドブックの重複データの排除の章を参照してください）。同意（インフォームド・コンセント）が取れていない記録（“concent” 変数（列）が false）も除外する必要があります。

下記例では、除外したいケースを抽出して、別のデータフレーム（dropped）に保存します。この抽出操作により、調査分析から除外された症例を明示できます。次に、dplyr パッケージの anti_join() 関数を使用して、調査データ（survey_data）オブジェクトからこれらの除外された症例を削除します。

危険： 体重の変数や、調査デザインに関連する変数（年齢、性別、層やクラスターの変数など）に欠損値があってはいけません。

## 除外した症例を別のデータフレームとして保存しておけば、除外理由を説明できる（例：同意がないもしくは、誤った市区町村・クラスター)
dropped <- survey_data %>% 
  filter(!consent | is.na(startdate) | is.na(enddate) | village_name == "other")

## 除外された症例データセットを使用して、調査データセットから未使用の記録（行）を削除する
survey_data <- anti_join(survey_data, dropped, by = names(dropped))

前述のように、3 つの異なる標本抽出法（層化抽出法、クラスター抽出法、層化クラスター抽出法）について、重みを追加する方法を示します。これらの抽出法には、母集団、標本の両方とも，またはいずれか一方の集団に関する情報が必要です。下記例では、層化クラスター抽出法のための重みを追加するコードを使用しますが、研究デザインに最も適したものを使用してください。

# 層化抽出法 ------------------------------------------------------------------------
# 年齢層別（age_group）、性別（sex）、担当保健所（health_district）各変数（列）についての各個人の重みを含んでいる（層化されている）
# "surv_weight_strata" という変数（列）を作成する
survey_data <- add_weights_strata(x = survey_data,
                                         p = population,
                                         surv_weight = "surv_weight_strata",
                                         surv_weight_ID = "surv_weight_ID_strata",
                                         age_group, sex, health_district)

## クラスター抽出法 ------------------------------------------------------------------

# 1 世帯あたりのインタビュー対象者（interviewed）の人数を得る
# 次に、世帯（親クラスター）の識別子を格納する変数のカウントを持つ変数（列）を追加する（cluster_counts）
survey_data <- survey_data %>%
  add_count(index, name = "interviewed")


## クラスターの重みを作成する
survey_data <- add_weights_cluster(x = survey_data,
                                          cl = cluster_counts,
                                          eligible = member_number,
                                          interviewed = interviewed,
                                          cluster_x = village_name,
                                          cluster_cl = cluster,
                                          household_x = index,
                                          household_cl = households,
                                          surv_weight = "surv_weight_cluster",
                                          surv_weight_ID = "surv_weight_ID_cluster",
                                          ignore_cluster = FALSE,
                                          ignore_household = FALSE)


# 層化クラスター抽出法 --------------------------------------------------------------
# クラスターと層化をかけ合わせた重み（surv_weight_strata * surv_weight_cluster）を作成する（surv_weight_cluster_strata）
survey_data <- survey_data %>%
  mutate(surv_weight_cluster_strata = surv_weight_strata * surv_weight_cluster)

26.6 調査デザインオブジェクト

研究デザイン（抽出法）に合わせて調査デザインオブジェクト（survey design object）を作成します。データフレーム（survey_data）と同様に、重み付け人口比率などの計算に使用します。必要な変数（列）がすべて作成されていることを確認してから行います。

抽出法の設定には下記の 4 つの選択肢がありますが、使用しない抽出法に関するコードはコメントアウトしてください。 - 単純無作為抽出法 - 層化抽出法 - クラスター抽出法 - 層化クラスター抽出法

このテンプレートでは、2 つの異なる層（担当保健所 A と B）で標本調査をクラスター化したと仮定します。したがって、全体の推定値を得るためには、クラスターと層の重みを組み合わせる必要があります。

前述したように、抽出法に沿った分析を行うために 2 つのパッケージが用意されています。古典的なものは survey パッケージで、それから tidyverse コーディングに適したオブジェクトや関数を作る srvyr パッケージというラッパーがあります。両方ともデモを行いますが、本章のコードのほとんどは srvyr パッケージベースのオブジェクトを使用することに注意してください。1 つの例外は、gtsummary パッケージが survey パッケージのオブジェクトしか受け付けないことです。

26.6.1 Survey パッケージ

survey パッケージの構文は、base R の構文を効果的に使用するため、パイプ（%>%）やその他の dplyr パッケージの構文を使用することはできません。survey パッケージでは、svydesign() 関数を使用して、適切なクラスタ、重み、層を持つ調査オブジェクト（survey object）を定義します。

注釈：svydesign() の引数に指定する変数の前にチルダ記号（~）を記述しなければいけません。なぜなら、survey パッケージは、formula オブジェクトに基づいて関数の引数に変数を割り当てる base R の構文を使用しているためです。

# 単純無作為抽出法 ------------------------------------------------------------------
base_survey_design_simple <- svydesign(ids = ~1, # クラスター ID なしを意味する 1 を指定
                   weights = NULL,               # 重みはなし
                   strata = NULL,                # 単純無作為抽出法であるため層の指定はなし
                   data = survey_data            # データセットを指定する必要あり
                  )

## 層化抽出法 -----------------------------------------------------------------------
base_survey_design_strata <- svydesign(ids = ~1,  # クラスター ID なしを意味する 1 を指定
                   weights = ~surv_weight_strata, # 以前のコードチャンクで作成した層化に対応する重みの変数を指定
                   strata = ~health_district,     # 層化抽出法の指定は担当保健所（health_district）ごとに層化する
                   data = survey_data             # データセットを指定する必要あり
                  )

# クラスター抽出法 ------------------------------------------------------------------
base_survey_design_cluster <- svydesign(ids = ~village_name, # クラスター ID を市区町村（village_name）に指定
                   weights = ~surv_weight_cluster, # 以前のコードチャンクで作成したクラスターに対応する重みの変数を指定
                   strata = NULL,                 # クラスター抽出法であるため層の指定はなし
                   data = survey_data              # データセットを指定する必要あり
                  )

# 層化クラスター抽出法 --------------------------------------------------------------
base_survey_design <- svydesign(ids = ~village_name,      # クラスター ID に市区町村（village_name）を指定
                   weights = ~surv_weight_cluster_strata, # 以前のコードチャンクで作成した重みの変数を指定
                   strata = ~health_district,             # 層化抽出法の指定は担当保健所（health_district）で層化を行う
                   data = survey_data                     # データセットを指定する必要あり
                  )

26.6.2 Srvyr パッケージ

Srvyr パッケージでは、as_survey_design() 関数を使うことができます。as_survey_design() は、svydesign() と同じ引数を取りますが、パイプ演算子（%>%）を使えるので、チルダ記号（~）を使う必要はありません。

## 単純無作為抽出法 ------------------------------------------------------------------
survey_design_simple <- survey_data %>% 
  as_survey_design(ids = 1, # クラスター ID なしを意味するは 1 を指定 
                   weights = NULL, # 重みはなし
                   strata = NULL # 単純無作為抽出法であるため層の指定はなし
                  )
## 層化抽出法 ------------------------------------------------------------------------
survey_design_strata <- survey_data %>%
  as_survey_design(ids = 1, # クラスター ID なしは 1 を指定
                   weights = surv_weight_strata, # 以前のコードチャンクで作成した重みの変数を指定
                   strata = health_district # 層化抽出法の指定は担当保健所（health_district）ごとに層化する
                  )
## クラスター抽出法 ------------------------------------------------------------------
survey_design_cluster <- survey_data %>%
  as_survey_design(ids = village_name, # クラスター ID を指定
                   weights = surv_weight_cluster, # 以前のコードチャンクで作成した重みの変数を指定
                   strata = NULL # 抽出は単純に（層はなし）
                  )

## 層化クラスター --------------------------------------------------------------
survey_design <- survey_data %>%
  as_survey_design(ids = village_name, # クラスター ID に市区町村（village_name）を指定
                   weights = surv_weight_cluster_strata, # 以前のコードチャンクで作成した重みの変数を指定
                   strata = health_district # 層化抽出法の指定は担当保健所（health_district）ごとに層化する
                  )

26.7 記述的分析

基本的な記述的分析と視覚化は、ハンドブックの他の章で広範囲にカバーされているので、ここでは触れません。詳細は、記述統計表の作り方、簡単な統計的検定、見やすい表の作り方、ggplot の基礎、R Markdown で作るレポートの章を参照してください。

このセクションでは、標本の偏りを調査し、その偏りの程度を視覚化する方法に焦点を当てます。また、沖積図・サンキーダイアグラムを使って、調査環境における人口の流れを視覚化することも検討します。

一般的には、以下のような記述的な分析を含めることを検討すべきです。

対象となるクラスター、世帯、個人の最終的な数
分析から除外された人の数とその理由
クラスターあたりの世帯数と世帯あたりの人数の中央値（範囲）

26.7.1 標本抽出のバイアス

各年齢層の割合を、標本と母集団の間で比較します。これは、潜在的な標本抽出のバイアスを明らかにするために重要です。同様の操作を性別に対して適用して分布を調べられます。

標本と母集団の比較（二項検定）で設定される p 値は単なる指標であり、母集団と比較した標本集団の分布の記述的な議論（または次のセクションの年齢ピラミッドによる視覚化）は、二項検定自体よりも重要であることに注意してください。なぜなら、標本サイズを大きくすると、データを重み付けした後では、p 値は無関係な差になることが多いからです。

## 標本の人数と年齢層ごとの人口構成比率
ag <- survey_data %>% 
  group_by(age_group) %>% 
  drop_na(age_group) %>% 
  tally() %>% 
  mutate(proportion = n / sum(n), 
         n_total = sum(n))

## 母集団の人数と年齢層ごと人口構成比率
propcount <- population %>% 
  group_by(age_group) %>%
    tally(population) %>%
    mutate(proportion = n / sum(n))

## 2 つの表の列を結合し、年齢別にグループ化し、二項検定を行って、
## 標本の n/total が母集団の n/total と有意に異なるかどうかを調べる
  ## 次の行では、2 つのデータセットの列の最後に接尾文字を追加している
joined_table <- left_join(ag, propcount, by = "age_group", suffix = c("", "_pop")) %>%
  group_by(age_group) %>%
  ## broom::tidy(binom.test()) は、二項検定の結果からデータフレームを作成し、
  ## 変数（列） p.value, parameter, conf.low, conf.high, method, and alternative を追加する
  ## ここでは p.value のみを使用する
  ## 信頼区間を報告したい場合は、他の列を含めることができる
  mutate(binom = list(broom::tidy(binom.test(n, n_total, proportion_pop)))) %>%
  unnest(cols = c(binom)) %>% # important for expanding the binom.test data frame
  mutate(proportion_pop = proportion_pop * 100) %>%
  ## 偽陽性を補正するために p 値を調整する 
  ## （複数の年齢層をテストしたことによる）
  ## これは、多くの年齢カテゴリーがある場合にのみ違いが生じる
  mutate(p.value = p.adjust(p.value, method = "holm")) %>%
                      
  ## 0.001 を超える p 値のみ示す(0.001 未満は <0.001 として示す)
  mutate(p.value = ifelse(p.value < 0.001, 
                          "<0.001", 
                          as.character(round(p.value, 3)))) %>% 
  
  ## 列の名前を適切に変更する
  select(
    "Age group" = age_group,
    "Study population (n)" = n,
    "Study population (%)" = proportion,
    "Source population (n)" = n_pop,
    "Source population (%)" = proportion_pop,
    "P-value" = p.value
  )

26.7.2 人口統計ピラミッド

人口統計（または年齢-性別）ピラミッドは、母集団の分布を視覚化する簡単な方法です。また、調査の層別に年齢と性別の記述統計表を作成することも検討する価値があります。ここでは、上記で作成した調査デザインオブジェクト（study design object）を使用して、加重比例が可能な apyramid パッケージを使用して説明します。人口ピラミッドとリッカート尺度の章で、人口統計ピラミッドを作成するための他の選択肢について詳しく説明しています。また、apyramid パッケージの age_pyramid() ラッパー関数を使用して、人口比率のプロットを作成するためのコーディングを省きます。

前述の標本抽出のバイアスのセクションで例示した母集団と標本の人口比率の差の二項検定のように、このセクションでは、我々の標本抽出した集団が元の母集団と実質的に異なるかどうか、そして重み付けがこの差を修正するかどうかを可視化を扱います。コードを実行するために、ggplot パッケージの出力を並べて表示するために patchwork パッケージを使用します。詳細はハンドブックの ggplot のヒントの章の「図の中に図を挿入する」のセクションを参照してください。このセクションでは、「母集団」、「重み付けされていない標本集団」、「重み付けされた標本集団」を視覚化します。また、調査の各層ごとに視覚化も可能です。ここでの例では、stack_by = "health_district" という引数を使用します（詳細は ?plot_age_pyramid を参照してください）。

注釈：人口ピラミッドのプロットでは、x 軸と y 軸が反転します。

## x 軸の範囲とラベルを定義する ---------------------------------------------
## (これらの数値をグラフの値として更新する)
max_prop <- 35      # 表示させたい人口比率の上限を選ぶ 
step <- 5           # 表示させたいラベル間のスペースを選ぶ

## 上記の数値を使って、軸の区切りのベクトルを定義する
breaks <- c(
    seq(max_prop/100 * -1, 0 - step/100, step/100), 
    0, 
    seq(0 + step / 100, max_prop/100, step/100)
    )

## 上記の数値を使って、軸の範囲のベクトルを定義する
limits <- c(max_prop/100 * -1, max_prop/100)

## 上記の数値を使って、軸のラベルのベクトルを定義する
labels <-  c(
      seq(max_prop, step, -step), 
      0, 
      seq(step, max_prop, step)
    )


## プロットを個別に作成する  --------------------------------------------------

## 母集団をプロットする
## 注：全体の人口を対象にして集計する必要がある（例：担当保健所によるクラスタリングを行わない）
source_population <- population %>%
  ## 年齢と性別が因子型であることを明示する
  mutate(age_group = factor(age_group, 
                            levels = c("0-2", 
                                       "3-14", 
                                       "15-29",
                                       "30-44", 
                                       "45+")), 
         sex = factor(sex)) %>% 
  group_by(age_group, sex) %>% 
  ## 各担当保健所内の人口数を足し合わせる
  summarise(population = sum(population)) %>% 
  ## グループ化を解除して全体の比率を算出する
  ungroup() %>% 
  mutate(proportion = population / sum(population)) %>% 
  ## 人口ピラミッドをプロットする
  age_pyramid(
            age_group = age_group, 
            split_by = sex, 
            count = proportion, 
            proportional = TRUE) +
  ## x 軸（表示状は垂直方向に表示される）のラベルのみを表示（設定しないと、他の 3 つのプロットすべてにラベル表示が繰り返される）
  labs(title = "Source population", 
       y = "", 
       x = "Age group (years)") + 
  ## すべてのプロットの x 軸（垂直表示）のスケールを同一にする 
  scale_y_continuous(breaks = breaks, 
    limits = limits, 
    labels = labels)
  
  
## 重み付けされていない標本集団をプロットする
sample_population <- age_pyramid(survey_data, 
                 age_group = "age_group", 
                 split_by = "sex",
                 proportion = TRUE) + 
  ## y 軸（表示状は水平方向に表示される）のラベルのみを表示（設定しないと、他の 3 つのプロットすべてにラベル表示が繰り返される）
  labs(title = "Unweighted sample population", 
       y = "Proportion (%)", 
       x = "") + 
  ## すべてのプロットの y 軸（水平表示）のスケールを同一にする 
  scale_y_continuous(breaks = breaks, 
    limits = limits, 
    labels = labels)


## 重み付けされた標本人口をプロットする
weighted_population <- survey_design %>% 
  ## make sure the variables are factors
  mutate(age_group = factor(age_group), 
         sex = factor(sex)) %>%
  age_pyramid(
    age_group = "age_group",
    split_by = "sex", 
    proportion = TRUE) +
  ## グラフタイトルを表示（設定しないと、他の 3 つのプロットすべてにラベル表示が繰り返される）
  labs(title = "Weighted sample population", 
       y = "", 
       x = "")  + 
  ## すべてのプロットの y 軸（水平表示）のスケールを同一にする 
  scale_y_continuous(breaks = breaks, 
    limits = limits, 
    labels = labels)

## 3 つのプロットを組み合わせる  ----------------------------------------------------
## 隣接する 3 つのプロットを + で結合する
source_population + sample_population + weighted_population + 
  ## 凡例を 1 つだけ表示し、テーマを定義する 
  ## テーマを plot_layout() と組み合わせるための & の使用法に注意する
  plot_layout(guides = "collect") & 
  theme(legend.position = "bottom",                    # move legend to bottom
        legend.title = element_blank(),                # remove title
        text = element_text(size = 18),                # change text size
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1) # turn x-axis text
       )

26.7.3 沖積図・サンキーダイアグラム

個人個人の観測開始とアウトカムを視覚化することは、概要を把握するのに非常に役立ちます。移動人口比率を考慮する非常に分かりやすいグラフ描画方法があります。また、コホートや、個人の状態に遷移があるその他の状況などを描画する数多くのグラフがあります。これらのグラフには、沖積図、サンキー、パラレルセットなど、いくつかの異なる表現があり、詳細はハンドブックのフローチャート・サンキー図・タイムラインの章に記載されています。

## 標本調査データの要約
flow_table <- survey_data %>%
  count(startcause, endcause, sex) %>%  # 観察期間開始理由（startcause）、観察期間終了理由（endcause）、性別（sex）変数（列）のカウントを取得する
  gather_set_data(x = c("startcause", "endcause"))     # プロットのために startcause、endcause の形式を縦長から横長に変更する


## データセットをプロットする
  ## x 軸は観察期間開始理由（startcause）と観察期間終了理由（endcause）
  ## gather_set_data では、可能な組み合わせごとに ID を生成する
  ## y で分割すると、可能性のある開始と終了の組み合わせを取得できる
  ## 値を n とすると、数となる（比率に変更することも可能）
ggplot(flow_table, aes(x, id = id, split = y, value = n)) +
  ## 性別での色分け 
  geom_parallel_sets(aes(fill = sex), alpha = 0.5, axis.width = 0.2) +
  ## ラベルボックスをグレーにする
  geom_parallel_sets_axes(axis.width = 0.15, fill = "grey80", color = "grey80") +
  ## 文字の色や角度の変更する（調整が必要）
  geom_parallel_sets_labels(color = "black", angle = 0, size = 5) +
  ## 軸ラベルを削除する
  theme_void()+
  # 凡例を図の下に移動
  theme(legend.position = "bottom")

26.8 重み付けされた人口比率

このセクションでは、重み付けされた人口数と人口比率の表を、関連する信頼区間とデザイン効果（標本の分散の、同じ要素数の単純無作為標本の分散に対する比）とともに作成する方法を詳しく説明します。survey パッケージ、srvyr パッケージ、sitrep パッケージ、gtsummary パッケージの関数を使った 4 つの異なる方法があります。標準的な疫学スタイルの表を作成するための最小限のコーディングには、srvyr パッケージのラッパーである sitrep 関数をお勧めしますが、これはまだ CRAN に掲載されておらず、将来変更される可能性があることに注意してください。それ以外では、srvyr パッケージが tidyverse のワークフローに最もうまく適合するのに対し、survey パッケージのコードは長期的に最も安定していると思われます。gtsummary パッケージの関数群は多くの可能性を秘めていますが、この記事を書いている時点では実験的で不完全なもののようです。

26.8.1 Survey パッケージ

survey パッケージから svyciprop() 関数を使用して、重み付けされた人口比率とそれに伴う 95% 信頼区間を得ることができます。適切なデザイン効果は、svyprop() ではなく、svymean() を使用して抽出できます。svyprop()は、0 から 1 の間の値（または TRUE/FALSE）しか受け付けませんので、カテゴリー変数は使えないことに注意する必要があります。

注釈： survey パッケージ由来の関数は、srvyr パッケージデザイン・オブジェクトも受け入れますが、ここでは一貫性を保つために survey パッケージデザイン・オブジェクトを使用しています。

## 重み付けされた死亡数を生成する
svytable(~died, base_survey_design)

## died
##      FALSE       TRUE 
## 1406244.43   76213.01

## 重み付けされた人口比率を生成する
svyciprop(~died, base_survey_design, na.rm = T)

##               2.5% 97.5%
## died 0.0514 0.0208  0.12

## デザイン効果を得る
svymean(~died, base_survey_design, na.rm = T, deff = T) %>% 
  deff()

## diedFALSE  diedTRUE 
##  3.755508  3.755508

上記のような survey パッケージの関数を組み合わせて、以下のように独自に定義した svy_prop という関数を作成できます。この関数と purrr パッケージの map() を使って、複数の変数を反復処理し、表を作成できます。purrr パッケージの詳細については、ハンドブックのループと反復処理・リストの操作の章を参照してください。

# 重み付けされた死亡数、人口比率、信頼区間、デザイン効果を計算する関数を定義する
# x は引用符で囲まれた変数
# design は survey パッケージのデザインオブジェクト

svy_prop <- function(design, x) {
  
  ## 興味のある変数を計算式に入れる 
  form <- as.formula(paste0( "~" , x))
  ## svytable からカウントの TRUE 列だけを残す
  weighted_counts <- svytable(form, design)[[2]]
  ## 割合を計算する（100 倍にして ％ を算出）
  weighted_props <- svyciprop(form, design, na.rm = TRUE) * 100
  ## 信頼区間を抽出し、乗算してパーセンテージを求める
  weighted_confint <- confint(weighted_props) * 100
  ## svymean を使ってデザイン効果を計算し、TRUE の列だけを残す
  design_eff <- deff(svymean(form, design, na.rm = TRUE, deff = TRUE))[[TRUE]]
  
  ## 1 つのデータフレームにまとめる
  full_table <- cbind(
    "Variable"        = x,
    "Count"           = weighted_counts,
    "Proportion"      = weighted_props,
    weighted_confint, 
    "Design effect"   = design_eff
    )
  
  ## テーブルをデータフレームとして返す
  full_table <- data.frame(full_table, 
             ## remove the variable names from rows (is a separate column now)
             row.names = NULL)
  
  ## 数値を数字型に戻す
  full_table[ , 2:6] <- as.numeric(full_table[, 2:6])
  
  ## データフレームを返す
  full_table
}

## 複数の変数を反復してテーブルを作成する 
purrr::map(
  ## 関心のある変数を定義する
  c("left", "died", "arrived"), 
  ## 使用する関数とその関数の引数（design）を指定する
  svy_prop, design = base_survey_design) %>% 
  ## リストを単一のデータフレームにまとめる
  bind_rows() %>% 
  ## 四捨五入する
  mutate(across(where(is.numeric), round, digits = 1))

##   Variable    Count Proportion X2.5. X97.5. Design.effect
## 1     left 701199.1       47.3  39.2   55.5           2.4
## 2     died  76213.0        5.1   2.1   12.1           3.8
## 3  arrived 761799.0       51.4  40.9   61.7           3.9

26.8.2 Srvyr パッケージ

srvyr パッケージでは、dplyr パッケージの構文を使って表を作成できます。survey_mean() が使用され、割合引数が指定されていること、また、デザイン効果の計算にも同じ関数が使用されていることに注意してください。これは、srvyr パッケージが、前述のセクションで使われている survey パッケージ関数群の svyciprop() と svymean() の両方を内包しているからです。

注釈： srvyr パッケージを使ってカテゴリカル変数から人口比率を求めることはできないようです。もしこの操作が必要であれば、sitrep パッケージを使用する以降のセクションをチェックしてください。

## srvyr パッケージのデザイン・オブジェクトの使用
survey_design %>% 
  summarise(
    ## 重み付けされた死亡数を生成する
    counts = survey_total(died), 
    ## 重み付けされた人口比率と信頼区間の生成する 
    ## 100 倍にしてパーセンテージを算出する
    props = survey_mean(died, 
                        proportion = TRUE, 
                        vartype = "ci") * 100, 
    ## デザイン効果を生成する
    deff = survey_mean(died, deff = TRUE)) %>% 
  ## 興味のある行だけを残す
  ## (標準誤差をドロップして、人口比率の計算を繰り返す)
  select(counts, props, props_low, props_upp, deff_deff)

## # A tibble: 1 × 5
##   counts props props_low props_upp deff_deff
##    <dbl> <dbl>     <dbl>     <dbl>     <dbl>
## 1 76213.  5.14      2.08      12.1      3.76

ここでも purrr パッケージを使って、複数の変数を反復処理する関数を書くことができます。purrr の詳細については、ハンドブックのループと反復処理・リストの操作の章を参照してください。

# 重み付けされた死亡数、人口比率、信頼区間、デザイン効果を計算する関数の定義する
# design は survey パッケージのデザインオブジェクト
# x は引用符で囲まれた変数


srvyr_prop <- function(design, x) {
  
  summarise(
    ## survey パッケージのデザインオブジェクトを使用する
    design, 
    ## 重み付けされた死亡数を生成する 
    counts = survey_total(.data[[x]]), 
    ## 重み付けされた人口比率と信頼区間の生成する
    ## 100 倍にしてパーセンテージを求める 
    props = survey_mean(.data[[x]], 
                        proportion = TRUE, 
                        vartype = "ci") * 100, 
    ## デザイン効果を求める
    deff = survey_mean(.data[[x]], deff = TRUE)) %>% 
  ## 変数名を追加する
  mutate(variable = x) %>% 
  ## 興味のある行だけを残す
  ## (標準誤差をドロップして、人口比率の計算を繰り返す)
  select(variable, counts, props, props_low, props_upp, deff_deff)
  
}
  

## 複数の変数を反復してテーブルを作成する 
purrr::map(
  ## 関心のある変数を定義する
  c("left", "died", "arrived"), 
  ## 使用する関数とその関数の引数（design）を指定する
  ~srvyr_prop(.x, design = survey_design)) %>% 
  ## リストを単一のデータフレームにまとめる
  bind_rows()

## # A tibble: 3 × 6
##   variable  counts props props_low props_upp deff_deff
##   <chr>      <dbl> <dbl>     <dbl>     <dbl>     <dbl>
## 1 left     701199. 47.3      39.2       55.5      2.38
## 2 died      76213.  5.14      2.08      12.1      3.76
## 3 arrived  761799. 51.4      40.9       61.7      3.93

26.8.3 Sitrep パッケージ

sitrep パッケージの tab_survey() は、srvyr パッケージのラッパーであり、最小限のコーディングで重み付けされた表を作成できます。また、カテゴリカル変数の加重割合を計算することもできます。

## survey_design オブジェクトを使用する
survey_design %>% 
  ## 興味のある変数の名前を引用符をつけずに渡す
  tab_survey(arrived, left, died, education_level,
             deff = TRUE,   # デザイン効果を求める
             pretty = TRUE  # 人口比率と95%信頼区間を統合する
             )

## Warning: removing 257 missing value(s) from `education_level`

## # A tibble: 9 × 5
##   variable        value            n  deff ci               
##   <chr>           <chr>        <dbl> <dbl> <chr>            
## 1 arrived         TRUE       761799.  3.93 51.4% (40.9-61.7)
## 2 arrived         FALSE      720658.  3.93 48.6% (38.3-59.1)
## 3 left            TRUE       701199.  2.38 47.3% (39.2-55.5)
## 4 left            FALSE      781258.  2.38 52.7% (44.5-60.8)
## 5 died            TRUE        76213.  3.76 5.1% (2.1-12.1)  
## 6 died            FALSE     1406244.  3.76 94.9% (87.9-97.9)
## 7 education_level higher     171644.  4.70 42.4% (26.9-59.7)
## 8 education_level primary    102609.  2.37 25.4% (16.2-37.3)
## 9 education_level secondary  130201.  6.68 32.2% (16.5-53.3)

26.8.4 Gtsummary パッケージ

gtsummary パッケージには、信頼区間やデザイン効果を追加するための組み込み関数がまだないようです。ここでは、信頼区間を追加するための関数を定義し、tbl_svysummary() を使って作成した gtsummary パッケージのテーブルに信頼区間を追加する方法を示します。

confidence_intervals <- function(data, variable, by, ...) {
  
  ## 信頼区間を抽出し、乗算してパーセンテージを求める
  props <- svyciprop(as.formula(paste0( "~" , variable)),
              data, na.rm = TRUE)
  
  ## 信頼区間を抽出する
  as.numeric(confint(props) * 100) %>% ## 数値化してパーセンテージを乗算する
    round(., digits = 1) %>%           ## 1 桁に四捨五入する
    c(.) %>%                           ## 行列から数値を抽出する
    paste0(., collapse = "-")          ## 1 つの文字列にまとめる
}

## survey パッケージのデザインオブジェクトを使用する
tbl_svysummary(base_survey_design, 
               include = c(arrived, left, died),   ## 含めたい変数を定義する
               statistic = list(everything() ~ c("{n} ({p}%)"))) %>% ## 興味のある統計量を定義する
  add_n() %>%  ## 重み付けされた総数を加える
  add_stat(fns = everything() ~ confidence_intervals) %>% ## 信頼区間を加える
  ## カラムのヘッダーを変更する
  modify_header(
    list(
      n ~ "**Weighted total (N)**",
      stat_0 ~ "**Weighted Count**",
      add_stat_1 ~ "**95%CI**"
    )
    )

Characteristic	Weighted total (N)	Weighted Count¹	95%CI
arrived	1,482,457	761,799 (51%)	40.9-61.7
left	1,482,457	701,199 (47%)	39.2-55.5
died	1,482,457	76,213 (5.1%)	2.1-12.1
¹ n (%)

26.9 重み付けされた比率

重み付けされた比率（例えば死亡率のような）については同様に、survey パッケージや srvyr パッケージを使うことができます。さらに、複数の変数を反復処理する関数（前のセクションのものと同様）を書くこともできます。上記のように gtsummary パッケージのための関数を作れますが、現在のところ相当の処理を行う関数はパッケージにありません。

26.9.1 Survey パッケージ

ratio <- svyratio(~died, 
         denominator = ~obstime, 
         design = base_survey_design)

ci <- confint(ratio)

cbind(
  ratio$ratio * 10000, 
  ci * 10000
)

##       obstime    2.5 %   97.5 %
## died 5.981922 1.194294 10.76955

26.9.2 Srvyr パッケージ

survey_design %>% 
  ## 観察期間を考慮した死亡率
  summarise(
    mortality = survey_ratio(
      as.numeric(died) * 10000, 
      obstime, 
      vartype = "ci")
    )

## # A tibble: 1 × 3
##   mortality mortality_low mortality_upp
##       <dbl>         <dbl>         <dbl>
## 1      5.98         0.349          11.6

26.10 参考資料