20 データの欠損
この章では以下の内容について説明します：
- 欠損の評価
- 欠損による行のフィルタリング
- 時系列ごとの欠損のプロット
- プロットにおける
NAの処理の仕方 - 欠損値の補完：MCAR, MAR, MNAR
20.1 準備
パッケージの読み込み
以下のコードを実行すると、分析に必要なパッケージが読み込まれます。このハンドブックでは、パッケージを読み込むために、pacman パッケージの p_load() を主に使用しています。p_load() は、必要に応じてパッケージをインストールし、現在の R セッションで使用するためにパッケージを読み込む関数です。また、すでにインストールされたパッケージは、R の基本パッケージである base （以下、base R）の library() を使用して読み込むこともできます。R のパッケージに関する詳細は R basics の章をご覧ください。
pacman::p_load(
rio, # インポート／エクスポート
tidyverse, # データの管理と視覚化
naniar, # 欠損の評価と視覚化
mice # 欠損値の補完
)データのインポート
エボラ出血熱の流行をシミュレートしたデータセットをインポートします。お手元の環境でこの章の内容を実行したい方は、 クリックして「前処理された」ラインリスト（linelist）データをダウンロードしてください>（.rds 形式で取得できます）。データは rio パッケージの import() を利用してインポートしましょう（rio パッケージは、.xlsx、.csv、.rds など様々な種類のファイルを取り扱うことができます。詳細は、インポートとエクスポート の章をご覧ください。）
# ラインリストをインポートする
linelist <- import("linelist_cleaned.rds")最初の 50 行を以下に示します。
インポート時における欠損の変換
データをインポートする際に、欠損、と分類されるべき値については少し注意が必要です。例えば、99、999、“Missing”、空欄（““）、または 空のスペースを含むセル（” “） などです。インポートのコマンドと一緒にこれらの値を NA （R において欠損値を示す）に変換することが出来ます。
使用するコマンドがファイルのタイプによって異なるので、詳細については データの欠損 のインポートのセクションを参考にしてください。
20.2 R における欠損値
以下では、R における欠損の表示や評価の仕方、また、それに関連するいくつかの値や関数について見ていきたいと思います。
NA
R においては、欠損値は予約語（特別な値）として指定された NA で表示されます。これは引用符なしで用いることに注意が必要です。“NA” とは異なり、これは通常の文字と同じ扱いとなります（ビートルズの Hey Jude の歌詞と同じです）。
データ内では “99”、“Missing”、または “Unknown” というように欠損が表されていることもあります。また、空の文字値 ““（空欄のように見える）や、スペース” ” で表されていることもあります。これらについては、インポート時に NA に変換する もしくはデータを na_if() を用いてクリーニングする必要があります。
データクリーニングにおいて、逆に全ての NA を “Missing” に変換したり、replace_na() を用いたり、因子型の場合には fct_explicit_na() を用いたりして処理することもできます。
いろいろな NA
多くの場合は、NA が欠損値を示し、それで全てがうまくいきます。しかし、オブジェクトのデータ型（文字値、数値、など）によって異なる NA を用いる必要が出てくる場合があります。あまり頻繁に起こることではないですが、注意が必要です。
多いのは、dplyr パッケージの case_when() を用いて新しい列を作成する場合などです。 データクリーニングと主要関数 の章で述べたように、この関数はデータフレーム内の全ての行に対して、ある基準（コードの右辺に書かれている条件）に合致するかどうかを評価し、新しい値（コードの左辺に書かれている値）を作成します。右辺に書かれている値は、全て同じデータ型でないといけません。
linelist <- linelist %>%
# "age" 列から新しく "age_years" 列を作成する
mutate(age_years = case_when(
age_unit == "years" ~ age, # age が年であれば、そのままの値を返す
age_unit == "months" ~ age/12, # age が月であれば、12 で割る
is.na(age_unit) ~ age, # age のユニットが欠損であれば、years とする
TRUE ~ NA_real_)) # 他のすべての場合には欠損とする右辺に NA を持ってきたい場合は、以下に示している特殊な NA のうちのどれかを用いる必要があります。他の右辺の値が文字型の場合は “Missing” か、NA_character_ を使用する必要があります。他の右辺が日付型や数字型の場合は NA_real_ を用いります。条件式であった場合 NA を用いることが出来ます。
-
NA- 条件式 TRUE/FALSE において用いる -
NA_character_- 文字値に用いる -
NA_real_- 日付型、数字値に用いる
繰り返しになりますが、新しい行を case_when() を用いて作成しようとしない限りは、こうした例に出会うことはありません。より詳細について知りたい方は、R documentation on NA を参考にしてください。
NULL
他に R で用いられる値に、NULL があります。これは、条件が正でも誤でもないことを示しています。すなわち、値が定義されない関数などを適用した際に返されます。一般的には、関数を書いたり、ある特定の条件において NULL を返すことを指示する shiny app を書いたりする場合を除いて、NULL を値として指定することはありません。
Nullかどうかは is.null() を用いて評価することができ、as.null() を用いて変換することができます。
NULL と NA の違いについては、ブログのポスト を見てください。
NaN
NaN は値としてとることが不可能なものを表す特別な値です。例えば、0 を 0 で割れ、と R に指示したような場合です。is.nan() を用いて評価することが出来、また、補足関数としてis.infinite() や is.finite() が存在します。
Inf
Inf は無限大を表します。例えば、ある値を0で割った場合などです。
これらの欠損値がどのように影響するかを簡単な例で示しましょう。例えば、z <- c(1, 22, NA, Inf, NaN, 5) で作成される z というベクターまたはコラムがあるとします。
一番大きな値を見つけるために max() をこのコラムに適用することを考えます。na.rm = TRUE を用いることで、計算から NA を除くことが出来ます。しかし、Inf や NaN は残り、結果として Inf が返されます。これは、[ ] と is.finite() を用いて有限な値のみを含むサブセットを作成して計算することで解決することができます： max(z[is.finite(z)]) 。
例
| R コマンド | 出力 |
|---|---|
5 / 0 |
Inf |
0 / 0 |
NaN |
5 / NA |
NA |
5 / Inf |0NA - 5|NAInf / 5|Infclass(NA)| "logical"class(NaN)| "numeric"class(Inf)| "numeric"class(NULL)` |
“NULL” |
“NAs introduced by coercion” はよくある警告メッセージです。数値のベクターに文字値を挿入しようとするなどの不正な操作を行った場合に起こります。
as.numeric(c("10", "20", "thirty", "40"))## Warning: NAs introduced by coercion## [1] 10 20 NA 40NULL はベクターにおいては無視されます。
my_vector <- c(25, NA, 10, NULL) # define
my_vector # print## [1] 25 NA 10ある一つの値の分散を求めようとすると NA が返されます。
var(22)## [1] NA20.3 便利な関数
以下では、R の base R の関数で、欠損値を評価したり扱ったりする際に便利なものを示します。
is.na() と !is.na()
is.na() は欠損値を特定する場合、または欠損でない値を特定する場合（! を一緒に用いります）に使います。どちらも理論値（TRUE or FALSE）が返されます。 返されたベクトルに sum() を適用、例えば sum(is.na(linelist$date_outcome)) とすることで、TRUE の数を計算することが出来ます。
## [1] FALSE FALSE FALSE TRUE FALSE TRUE FALSE
!is.na(my_vector)## [1] TRUE TRUE TRUE FALSE TRUE FALSE TRUE## [1] 2
na.omit()
この関数をデータフレームに対して適用することで、一つでも欠損値を含む行を除くことが出来ます。これもまた、base R に含まれます。
ベクトルに対して適用すると、NA を除くことが出来ます。例えば：
na.omit(my_vector)## [1] 1 4 56 5 22
## attr(,"na.action")
## [1] 4 6
## attr(,"class")
## [1] "omit"
drop_na()
データクリーニング の際に便利な tidyr パッケージの関数です。括弧内を空にして走らせると、欠損値を一つでも含む行を除くことが出来ます。括弧内で列の名前を指定することで、その列内の値が欠損している行を除くことが出来ます。列の指定には “tidyselect” を用いることもできます。
20.4 データフレーム内の欠損を評価する
データフレーム linelist 内の欠損を評価し、視覚化するのには naniar パッケージが便利です。
# パッケージをインストールするもしくは読み込む
pacman::p_load(naniar)欠損の数を定量化する
全ての欠損の値の割合を示すには pct_miss() を、数を示すには n_miss() を用いります。
# すべてのデータフレーム内の欠損値の割合
pct_miss(linelist)## [1] 6.688745以下の 2 つの関数は、それぞれ、欠損値を含む行の割合、または欠損が一つもない行の割合を返してくれます。NA は欠損を意味しますが、"" や " " は欠損としてはカウントされないことに注意してください。
# 欠損値を含む行の割合を示す
pct_miss_case(linelist) # 数を示すには n_miss() を用いる## [1] 69.12364
# コンプリートケース（欠損のない行）の割合を示す
pct_complete_case(linelist) # 数を示すには n_complete() を用いる## [1] 30.87636欠損状況を視覚化する
gg_miss_var() を用いることで、それぞれの列における欠損値の数（もしくは割合）を示すことが出来ます。
- プロットをグループごとに見たい場合は、
facet =で列の名前（引用符にいれない）を指定することが出来ます - デフォルトでは、割合ではなく数が表示されます。割合にしたい場合は
show_pct = TRUEで指定します - 通常の
ggplot()通り、軸やタイトルのラベルを追加したい場合は+ labs(...)を用いて行うことが出来ます
gg_miss_var(linelist, show_pct = TRUE)
以下では、データが %>% を用いて関数に渡されています。facet = を用いても、データを分けることが出来ます。
linelist %>%
gg_miss_var(show_pct = TRUE, facet = outcome)
どの値が欠損しているかを示すために、vis_miss() を用いてデータフレームをヒートマップで表すことが出来ます。また 特定の列をデータフレームから select() して関数にその列だけを渡すこともできます。
# データフレーム全体に対して欠損のヒートプロットを示す
vis_miss(linelist)
欠損状況の関係性を探り、可視化する
存在しないものをどうやって視覚化すればよいのでしょうか？デフォルトでは、ggplot() を使うと、欠損値を除いてプロットを作成してしまいます。
naniar パッケージの geom_miss_point() を用いれば解決します。2 つの列から散布図を作成するときに、ある一方の変数は存在し、片方の変数は欠損である場合には、欠損値をその列の最小値よりもさらに 10% 小さな値に変換し、違う色で散布図上に示してくれます。
以下に示している散布図では、赤いドットはある一方の変数は存在しているが片方の変数が欠損している場合を示しています。これによって欠損していない値に対する欠損している値の分布を視覚的に確認することが出来ます。
別の列によって層別化したあとに欠損を評価したい場合は、gg_miss_fct() を使います。これは、因子型（または日付）の列の値ごとにデータフレームの欠損の割合をヒートマップで示してくれます。
gg_miss_fct(linelist, age_cat5)
この関数を日付ごとにデータフレームを示すのに用いると、欠損が時系列を追うごとにどのように変化しているかを示すことも出来ます。
gg_miss_fct(linelist, date_onset)## Warning: Removed 29 rows containing missing values (`geom_tile()`).
「付随する」列
ある一方の列の欠損を、もう片方の列の値ごとに示す別の方法として、naniar パッケージの作成する「付随する列」を用いるものがあります。bind_shadow() を用いると、全ての列に対して NA または not NA の 2 値変数を作成し、“_NA” を末尾に付けた新しい列として元のデータセットに結合してくれます。すなわち、データセットの列は 2 倍になります。
## [1] "case_id" "generation"
## [3] "date_infection" "date_onset"
## [5] "date_hospitalisation" "date_outcome"
## [7] "outcome" "gender"
## [9] "age" "age_unit"
## [11] "age_years" "age_cat"
## [13] "age_cat5" "hospital"
## [15] "lon" "lat"
## [17] "infector" "source"
## [19] "wt_kg" "ht_cm"
## [21] "ct_blood" "fever"
## [23] "chills" "cough"
## [25] "aches" "vomit"
## [27] "temp" "time_admission"
## [29] "bmi" "days_onset_hosp"
## [31] "case_id_NA" "generation_NA"
## [33] "date_infection_NA" "date_onset_NA"
## [35] "date_hospitalisation_NA" "date_outcome_NA"
## [37] "outcome_NA" "gender_NA"
## [39] "age_NA" "age_unit_NA"
## [41] "age_years_NA" "age_cat_NA"
## [43] "age_cat5_NA" "hospital_NA"
## [45] "lon_NA" "lat_NA"
## [47] "infector_NA" "source_NA"
## [49] "wt_kg_NA" "ht_cm_NA"
## [51] "ct_blood_NA" "fever_NA"
## [53] "chills_NA" "cough_NA"
## [55] "aches_NA" "vomit_NA"
## [57] "temp_NA" "time_admission_NA"
## [59] "bmi_NA" "days_onset_hosp_NA"これらの「付随する列」は、他のどの列に対してでも、欠損している値の割合をプロットするのに使えます。
例えば、以下のプロットでは、date_hospitalisation 列の値ごとに、days_onset_hosp 列（発症してから入院までの日数）における欠損の割合を示しています。x 軸の列の密度をプロットしており、興味のある付随する列によって層別化（color =）していることに注意してください。x 軸が数値もしくは日付の列である場合、このプロットはうまく実行されます。
ggplot(data = shadowed_linelist, # 付随する列を含むデータフレーム
mapping = aes(x = date_hospitalisation, # 数値もしくは日付の列
colour = age_years_NA)) + # 興味のある付随する列
geom_density() # 密度曲線をプロットする
以下に示すように、統計量のサマリを層別化して示すのにも「付随する列」を使用することが出来ます。
linelist %>%
bind_shadow() %>% # 「付随する列」を作成する
group_by(date_outcome_NA) %>% # 層別化に用いる「付随する列」を指定する
summarise(across(
.cols = age_years, # 統計量を示したい変数を指定する
.fns = list("mean" = mean, # 興味のある統計量を指定する
"sd" = sd,
"var" = var,
"min" = min,
"max" = max),
na.rm = TRUE)) # 統計量の計算に用いるその他のコマンド## # A tibble: 2 × 6
## date_outcome_NA age_years_mean age_years_sd age_years_var age_years_min
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 !NA 16.0 12.6 158. 0
## 2 NA 16.2 12.9 167. 0
## # ℹ 1 more variable: age_years_max <dbl>以下に、列の欠損値の割合を時系列ごとに示す他の方法を示します。naniar パッケージは使用しません。週ごとの観測値の欠損の割合を示しています。
- データを、興味のある単位時間（日、週、など）でまとめ、
NA（そしてその他興味のある値）を含む観測値の割合を計算します -
ggplot()を用いて、欠損の割合を折れ線グラフとしてプロットします
以下では、ラインリストを用いて、まず週ごとの値を新しい列として加え、値が欠損している週ごとの記録の割合を計算しています（7 日間での割合を計算したい場合は、以下のスクリプトとは少し異なります）。
outcome_missing <- linelist %>%
mutate(week = lubridate::floor_date(date_onset, "week")) %>% # 新しく week の列を作成する
group_by(week) %>% # 行を、week ごとにグルーピングする
summarise( # それぞれの値を week ごとにまとめる
n_obs = n(), # 観測数
outcome_missing = sum(is.na(outcome) | outcome == ""), # 欠損を含む観測数
outcome_p_miss = outcome_missing / n_obs, # 欠損を含む観測の割合
outcome_dead = sum(outcome == "Death", na.rm=T), # 死亡例の観測数
outcome_p_dead = outcome_dead / n_obs) %>% # 死亡例の観測の割合
tidyr::pivot_longer(-week, names_to = "statistic") %>% # ggplot のために week を除くすべての列をロング形式にピボットする
filter(stringr::str_detect(statistic, "_p_")) # 割合を示す値だけを残すこのデータを用いて、週ごとに欠損の割合を折れ線で示します。ggplot2 パッケージについてあまり詳しくない場合は、ggplotの基礎 の章を参考にしてください。
ggplot(data = outcome_missing)+
geom_line(
mapping = aes(x = week, y = value, group = statistic, color = statistic),
size = 2,
stat = "identity")+
labs(title = "Weekly outcomes",
x = "Week",
y = "Proportion of weekly records") +
scale_color_discrete(
name = "",
labels = c("Died", "Missing outcome"))+
scale_y_continuous(breaks = c(seq(0,1,0.1)))+
theme_minimal()+
theme(legend.position = "bottom")
20.5 欠損値を含むデータを使う
欠損値を含む行を取り除く
欠損値を含む行を除く簡単な方法は、dplyr パッケージの drop_na() を用いる方法です。
元の linelist には nrow(linelist) 行含まれています。欠損値を含む行を除いた後の行数は以下になります。
## [1] 1818また、ある特定の列内の欠損値を含む行を除くこともできます。
## [1] 5632順番に行を表示することも出来ますし、 “tidyselect” helper functions を使うこともできます。
## [1] 3029
ggplot() での NA の扱い方
たいていの場合、プロットから除かれた値の数はキャプションとして報告するのが賢明です。以下に例を示します。
ggplot() においては、labs() を追加し、その中で caption = としてキャプションを追加することが出来ます。キャプションの中で stringr パッケージの str_glue() を使うことで、文中に直接値を張り付けることが出来るので、データそれぞれの値を用いることが出来ます。
-
\nは、行替えに用いることが出来ます。 - もし、複数の列がプロットに表示されない値を含む場合（例えば、年齢あるいは性別がプロットに反映されている場合）、表示されていない数を正確に計算するためにはそれらの列で絞り込みを行う必要があります。
labs(
title = "",
y = "",
x = "",
caption = stringr::str_glue(
"n = {nrow(central_data)} from Central Hospital;
{nrow(central_data %>% filter(is.na(date_onset)))} cases missing date of onset and not shown.")) ggplot() コマンドを使用する前に、文字をオブジェクトとして保存し、str_glue() 内ではそのオブジェクト名を指定するだけの方が簡単な場合もあります。
因子型におけるNA
興味のある列が因子型である場合、NA を文字値にするには forcats パッケージの fct_explicit_na() を用いる必要があります。詳細については因子（ファクタ）型データ (#factors)の章を参考にしてください。デフォルトでは、新しい値は “(Missing)” として扱われますが、na_level = のコマンドによって変更することが出来ます。
pacman::p_load(forcats) # パッケージの読み込み
linelist <- linelist %>%
mutate(gender = fct_explicit_na(gender, na_level = "Missing"))
levels(linelist$gender)## [1] "f" "m" "Missing"20.6 欠損値の補完
単純に全ての欠損値を除いたデータセットを解析するのではなく、「ギャップを埋めて」、欠損値を補完してからデータを解析する方が良い場合もあります。単純にすべての欠損値を除くと、いろいろ問題があることがあります。例えば：
欠損値を含む観測値や大量の欠損値を含む変数を全て除くと、ある種の解析を行うパワーが減る可能性があります。例えば、先に示したように、linelist データにおいては、全ての変数のうち欠損値をまったく含まない観測は限られた割合しかありません。データセットの大部分を取り除いてしまうと、大量の情報を失ってしまうことになります！さらに、たいていの変数はそれなりの量の欠損値を含んでいます。そのため、ほとんどの分析において欠損値がたくさん含まれる変数を除去することはあまり適切ではありません。
欠損がなぜ生じているのかによっては、単に欠損を含まないデータのみで解析することは、解析結果にバイアスが生じたり、ミスリーディングな結果となることがあります。例えば、ある患者さんが熱や咳などの重要な症状の有無に関するデータを欠損しているとします。たとえば、明らかにそこまで重症でない患者さんについては熱や咳などの症状について記録をしていない可能性もあります。その場合は、単に欠損を含む観測値を除くと、データセット内の最も健康な人々を除くことになり、非常に結果にバイアスを生じることになります。
どのくらい欠損しているか、ということに加えて、なぜ欠損が生じてるのかを考えることが重要です。そうすることで、欠損データを補完することがどのくらい重要なのか、そしてどの方法で補完することが一番良いのかについて判断を下しやすくなるでしょう。
欠損データの種類
欠損データには一般的に、以下に挙げる 3 つの種類があります；
Missing Completely at Random (MCAR) これは、データが欠損する確率と、他のデータ内との変数との間に関係性がないことを意味します。すべてのケースで、欠損する確率がおなじこととなります。これは珍しい状況です。しかし、もしデータが MCAR であると信じる確固とした理由があるのであれば、欠損していないデータだけを補間を行わずに分析することは、結果にバイアスを生じさせません（とはいえ、パワーはいくらか失われます）。[TODO: consider discussing statistical tests for MCAR]
Missing at Random (MAR) この名前は少しミスリーディングで、というのも、MAR は、生じている欠損のパターンを他の変数の情報から規則的に予測できるタイプの欠損データのことだからです。例えば、全てのケースにおいて、熱のデータを欠損しているのは、実は寒気と痛みのある人は熱があると推測されたために熱を測らなかったためだとします。そうだとすると、寒気と痛みのあった観測値においては、発熱もまた生じていると容易に予測することが出来、欠損値を補完するためにこの情報を使用することが出来ます。実際には、これはもう少しスペクトラム的な話で、悪寒と痛みの両方を訴えており発熱のデータを欠損している患者は発熱もしている可能性が高いが、必ずしもそうというわけではない、ということです。しかし、完全に予測できるわけではなくとも、予測可能ではあります。欠損データの種類としては、最も多いタイプのものです。
Missing not at Random (MNAR) Not Missing at Random (NMAR) と呼ばれることもあります。このタイプは、欠損データの割合が系統的でなく、他の情報を用いて予測可能でもなく、また、ランダムに欠損が生じているわけでもないタイプの欠損のデータのことです。すなわち、データは何らかの明らかになっていない原因により欠損しているか、あるいは情報を持っていない何らかの規則に基づいて欠損しているといえます。例えば、年齢のデータが欠損しているのは、非常に高齢の患者は自分の年齢を知らないか、自分の年齢を開示したくなかったからだとします。この場合、年齢の欠損は、年齢データ自身に関連しており（すなわちランダムに生じているわけではない）、従って他の情報に基づいて欠損を予測することはできません。MNAR は複雑で、対処法としては、欠損を補完しようとするよりも、もっとデータを集めること、そしてなぜ欠損が生じているのかについての情報を収集することが良いでしょう。
まとめると、MCAR データを補完することは比較的容易ですが、一方、MNAR データを補完することは不可能ではないにしても非常にチャレンジングです。多くのデータの補完のメソッドは、MAR を仮定しています。
便利なパッケージ
Mmisc、missForest（ランダムフォレストによる欠損値の補完）、そして mice（連鎖方程式による多重代入法）などのパッケージは、欠損値の補完に便利です。このセクションでは、いろいろなテクニックを学ぶことの出来る mice パッケージについてのみ解説します。mice の開発者によるオンラインブックには、より詳細が示されています（https://stefvanbuuren.name/fimd/）。
mice パッケージを読み込みます；
pacman::p_load(mice)平均値代入法
例えば、非常に単純な解析をしている、あるいは MCAR を仮定する強い理由がある時には、数値の欠損を単にその変数の平均値で代入することが出来る場合があります。例えば、我々のデータセット内の気温の欠損値は MCAR であるか単に平年値である、と仮定できる可能性があります。気温の欠損値をデータセットの値の平均値で補完した新しい変数を作成するコードを示します。しかし、多くの場合では単に平均値を補完するとバイアスを生じる可能性があるので、十分注意が必要です。
linelist <- linelist %>%
mutate(temp_replace_na_with_mean = replace_na(temp, mean(temp, na.rm = T)))同じような方法で、因子型の欠損値をある特定の値で補完することもできます。例えば、転帰について（“Death” または “Recover”）の情報が欠損している観測値はすべて死亡データであるとします（このデータセットについては、これは本当は正しくないのですが）。
linelist <- linelist %>%
mutate(outcome_replace_na_with_death = replace_na(outcome, "Death"))回帰代入法
いくらかアドバンスな方法として、何らかの統計モデルを用いて欠損値がどんな値らしいのかを推測し、それで補完する方法があります。ここでは、体温についての情報は欠損しているけれど年齢と発熱についての情報は存在する観測値に対して、年齢（年）と発熱を予測因子とした簡単な線形回帰モデルを用いて値を予測する方法を示しています。実際には、このような単純なモデルではなく、より正確なモデルを使用した方が良いでしょう。
simple_temperature_model_fit <- lm(temp ~ fever + age_years, data = linelist)
#単純な体温モデルを用いて、体温の欠損値を予測する
predictions_for_missing_temps <- predict(simple_temperature_model_fit,
newdata = linelist %>% filter(is.na(temp))) もしくは、体温の欠損値に対して同様にモデルを作成する方法として mice パッケージを用いる方法もあります。
model_dataset <- linelist %>%
select(temp, fever, age_years)
temp_imputed <- mice(model_dataset,
method = "norm.predict",
seed = 1,
m = 1,
print = F)## Warning: Number of logged events: 1
temp_imputed_values <- temp_imputed$imp$tempこれは、欠損値を予測値で補完する、という意味で missForest パッケージなどのよりアドバンスな方法と同じタイプの方法だと言えます。missForest の場合は、線形回帰ではなくランダムフォレストを用いて予測モデルを作成します。他のタイプのモデルを使用することもできます。この方法は、MCARの場合はうまくいきますが、MAR あるいは MNAR の方がデータをうまく説明する、と思っている場合には注意が必要です。補完の質は、予測モデルがどのくらい正確か、ということにも拠りますが、正確なモデルであっても補完されたデータのばらつきが少し過小となる場合もあります。
LOCF と BOCF
Last observation carried forward (LOCF) と baseline observation carried forward (BOCF) は、時系列および縦断データの補完方法です。考え方としては、欠損の生じている観測値よりも前に観測されたデータで欠損値を補完しようというものです。連続で値が欠損している場合は、最後に観測された値で補完します。
LOCF および BOCF のどちらも tidyr パッケージの fill() を用いて実装することが出来ます（HMISC、zoo、そして data.table などのパッケージを用いても行うことが出来ます）。fill() をどのように使うのか示すために、2000 年と 2001 年において四半期ごとの罹患数を含む単純な時系列データを作成します。しかし、Q1 以降の半期に関するデータは欠損しているので、補完する必要があります。データの縦横変換 の章においても、fill() の使用方法が示されています。
#単純なデータセットを作成する
disease <- tibble::tribble(
~quarter, ~year, ~cases,
"Q1", 2000, 66013,
"Q2", NA, 69182,
"Q3", NA, 53175,
"Q4", NA, 21001,
"Q1", 2001, 46036,
"Q2", NA, 58842,
"Q3", NA, 44568,
"Q4", NA, 50197)
#year の欠損値を補完する
disease %>% fill(year)## # A tibble: 8 × 3
## quarter year cases
## <chr> <dbl> <dbl>
## 1 Q1 2000 66013
## 2 Q2 2000 69182
## 3 Q3 2000 53175
## 4 Q4 2000 21001
## 5 Q1 2001 46036
## 6 Q2 2001 58842
## 7 Q3 2001 44568
## 8 Q4 2001 50197fill() を使用する前に、データが正しくソートされているかどうかを確認してください。fill() はデフォルトで「下向きに」補完しますが、 .direction パラメタを用いることで違う方向に補完することもできます。例えば、似たようなデータですが、今度は年の終わり（Q4）のみ記録があり、他は欠損しているデータセットがあるとします。
#先ほどと少しだけ違うデータセットを作成する
disease <- tibble::tribble(
~quarter, ~year, ~cases,
"Q1", NA, 66013,
"Q2", NA, 69182,
"Q3", NA, 53175,
"Q4", 2000, 21001,
"Q1", NA, 46036,
"Q2", NA, 58842,
"Q3", NA, 44568,
"Q4", 2001, 50197)
#year の欠損値を、今度は「上向き」に補完する
disease %>% fill(year, .direction = "up")## # A tibble: 8 × 3
## quarter year cases
## <chr> <dbl> <dbl>
## 1 Q1 2000 66013
## 2 Q2 2000 69182
## 3 Q3 2000 53175
## 4 Q4 2000 21001
## 5 Q1 2001 46036
## 6 Q2 2001 58842
## 7 Q3 2001 44568
## 8 Q4 2001 50197この例では、LOCF およびに BOCF は明らかに正しい方法でしたが、より複雑な状況ではこの方法が適切かどうかを判断するのが難しい場合もあります。例えば、初日から何日間か、患者の検査値が欠損しているとします。ある場合は、これは検査値に変動がないことを示していますが、患者が回復しており、値が初日とはまったく異なる可能性もあります！なので、少し注意してこれらのメソッドを使う必要があります。
多重代入法
先ほど紹介した mice パッケージの開発者によるオンラインブック（https://stefvanbuuren.name/fimd/）には多重代入法の説明と、なぜそれを使いたいのかが詳細に記されています。ここでは、方法の基本的な部分を説明します。
多重代入法では、尤もらしい値で欠損値を代入したデータセットを複数作成します（データによってはより多くもしくは少なくデータセットを作成したい場合もあると思いますが、mice パッケージのデフォルトでは、5 つのデータセットが作成されます）。違いとしては、代入される値は、ある単一の値ではなく、予測された分布に基づき特定される値（したがっていくらかのランダムさを含んでいます）である、という点です。それぞれのデータセットに対してある種の予測モデルを用いる（mice には様々なオプションがあり、例えば予測平均マッチング、ロジスティック回帰、ランダムフォレストなどがあります）点では同じですが、mice パッケージでは様々なモデルの詳細についてまで考慮することが出来ます。
代入されたデータセット群を作成したら、これらの新しいデータセットそれぞれに対して、もともと計画していた統計解析を行い、それらの結果を一つにプールします。この方法は、MCAR において、また MAR の多くの場合においてバイアスを減少するのに非常に有用で、より正確な標準誤差を得ることが出来ます。
以下では、多重代入法を用いて、ラインリストのデータセット（より単純にした model_dataset を用いります）で、年齢と発熱のデータから体温を予測する方法を示します。
# model_datasetにおける欠損をすべて代入し、新たに 10 つのデータセットを作成する
multiple_imputation = mice(
model_dataset,
seed = 1,
m = 10,
print = FALSE) ## Warning: Number of logged events: 1
model_fit <- with(multiple_imputation, lm(temp ~ age_years + fever))
base::summary(mice::pool(model_fit))## term estimate std.error statistic df p.value
## 1 (Intercept) 3.703143e+01 0.0270863456 1367.16240465 26.83673 1.583113e-66
## 2 age_years 3.867829e-05 0.0006090202 0.06350905 171.44363 9.494351e-01
## 3 feveryes 1.978044e+00 0.0193587115 102.17849544 176.51325 5.666771e-159ここでは、mice のデフォルトの予測モデルである予測平均マッチングで欠損を補完しました。そして、それぞれの代入されたデータセットを値の推定に用いり、線形回帰によりそれぞれのデータセットにおいて得られた推定値をプールします。ここでは触れませんでしたが、mice パッケージを用いる際にはより詳細な方法や多重代入法に関する様々な設定を考慮する必要があります。例えば、常に数値データであるとは限らず、その場合は他の代入方法を用いる必要があります（mice パッケージでは、他の様々なデータの種類や補完が出来ます）。欠損データが非常に問題となる場合は、より頑健な解析結果を得るためには、単にコンプリートケースによる解析を行うよりも、多重代入法を用いる方が良いと言えるでしょう。
20.7 参考資料
Vignette on the naniar package
Gallery of missing value visualizations
Online book about multiple imputation in R by the maintainer of the mice package