miceパッケージampute関数による欠測データの作成:(2)実行例

はじめに

今回は、以前の記事で概要を紹介したampute関数を用いて実際に欠測データを作成してみたいと思います。
amputeの概要については、よろしければ以前の記事もご覧ください。
mstour.hatenablog.com

使用するデータセット

まずは今回使用するパッケージとデータセットを読み込んでおきます。
ampute関数を使うためにmiceパッケージを、サンプルデータセットを得るためにlme4パッケージも読み込みます（後述）。
また、後ほどデータの加工やプロットに必要となるtidyverseも読み込んでおきます。
データセットは2種類を使ってみます。

library(tidyverse)
library(mice)
library(lme4)
set.seed(1923)

sampdata_1 = mice::boys
sampdata_2 = lme4::sleepstudy

1. boys データセット

1つめの例は、miceパッケージに入っている「boys」というデータセットです。
これはある年における少年の年齢・身長・体重などの特徴を調べた、cross-sectional dataと呼ばれる形式（大雑把に言うと、ある特定時点のスナップショット）のデータセットです。
このデータセットに対して、「年齢が高くなるほど欠測が起きやすい」という状況を想定して欠測を作りたいと思います。
ただし、boysデータセットにはすでに多くの欠測（NA）が含まれているので、このままではampute関数が使えません。また、文字型の変数もそのままでは使えないので、使いずらい変数や欠測レコードをあらかじめ削除しておくことにします。

sampdata_1r = sampdata_1 %>%
  select(-gen, -phb, -tv, -reg) %>%
  filter(complete.cases(.))

2行目で不要な変数を削除し、3行目でさらに欠測のないレコードだけを残す処理をしています。

2. sleepstudy データセット

2つめの例は、lme4パッケージに入っている「sleepstudy」データセットを用います。
これは睡眠時間制限を連日行い、その期間中に何らかの課題を実施した時の反応時間*1を毎日記録したもので、このようなデータセットはlongitudinal dataと呼ばれます。
この例では、「前日の測定値が当日欠測するか否かに影響する」というパターンを考えてみます。

なお、ampute関数では（というか、miceパッケージでは*2）データセットを横持ち（つまり、日々の測定値はそれぞれ別変数として保存され、横に広い形になる）にしておく必要があります。
sleepstudyデータセットは各個人の反復測定データが縦に並んだ（縦持ちの）形状なので、以下のように加工しておきます。

sampdata_2r = sampdata_2 %>%
  pivot_wider(names_from = Days, values_from = Reaction, names_prefix = "Day")

測定値の変数名として、測定日の数値の前に"Day"をつけることにしました。

事例1. boys データセット

早速1つめの例で欠測データを作っていきましょう。

前回ご紹介したように、ampute関数では欠測パターンや欠測割合など、いくつかのパラメータを指定する必要があります。

f:id:mstour:20210826071029j:plain

パラメータの与え方で少し手間取る可能性があるので、まずはいくつかのパラメータ情報の形式を確認してみましょう。
デフォルト設定でamputeを適用した後、設定を確認します。

defo_1 = ampute(sampdata_1r)
ptn_defo_1 = defo_1$patterns
frq_defo_1 = defo_1$freq
wgt_defo_1 = defo_1$weights

欠測パターンは、次のような行列型で定義されます。

> ptn_defo_1
  age hgt wgt bmi hc
1   0   1   1   1  1
2   1   0   1   1  1
3   1   1   0   1  1
4   1   1   1   0  1
5   1   1   1   1  0

各行が今回用いられた欠測パターンを表しており、0が欠測する変数、1が欠測しない変数を示します。
つまり、1行目のパターンは「ageのみが欠測する」、2行目は「hgtのみが欠測する」となり、合計5つの欠測パターンが存在していることになります。
自分で欠測パターンを定義するには、このルールにそって行列を作る必要があります。

各欠測パターンの起こる割合は、次のようなベクトル形式で表します。

> frq_defo_1
[1] 0.2 0.2 0.2 0.2 0.2

デフォルトでは、すべてのパターンの割合が等しくなるように設定されるようです。

各変数が欠測へ与える影響の大きさ（重み）は次のような行列形式になります。

> wgt_defo_1
  age hgt wgt bmi hc
1   0   1   1   1  1
2   1   0   1   1  1
3   1   1   0   1  1
4   1   1   1   0  1
5   1   1   1   1  0

この行列の1〜5行は、それぞれ欠測パターン1〜5に対応しています。
数値の大きさは欠測へ与える影響の大きさを示しており、0は欠測に全く関係しないことを意味しています。
したがって、1行目は欠測パターン1の「ageのみ欠測」にどの変数がどの程度影響するかを定義しており、デフォルトでは「age自身は欠測に影響せず、それ以外の変数は等しい影響度を持つ」という設定になっています*3。

そのほかのパラメータは基本的に数値をそのまま与えるだけなので、以上の点に注意すればよいでしょう。

では、ここからは設定を自分で決めて、欠測データを作ってみることにします。
最初に述べたように、年齢が高くなるほど欠測が多いようなデータを作ります。

まずは欠測パターンを定義します。
話を単純にするために、以下のような欠測パターンだけを考えることにします。

身長（hgt）だけが欠測する
身長（hgt）と体重（wgt）が欠測する

なお、本来はどちらの場合でもBMI（bmi）は計算できませんが、ここではそれは考えないことにします。
上記の2パターンの欠測を定義するには、次のように行列を用意すればOKです。

ptn_1 = matrix(c(1, 0, 1, 1, 1,
                 1, 0, 0, 1, 1),
               nrow = 2, byrow = T)
colnames(ptn_1) = colnames(sampdata_1r)

ptn_1
> ptn_1
     age hgt wgt bmi hc
[1,]   1   0   1   1  1
[2,]   1   0   0   1  1

次に各欠測パターンの割合は、身長のみが70%、身長・体重の両方が30%としておきます。

frq_1 = c(0.7, 0.3)

欠測に対する各変数の影響の重みについては、

年齢のみが欠測に影響する
欠測する変数自身の値は関係しない

という設定とします（欠測メカニズムで言うとMARとなります）。

wgt_1 = matrix(c(1, 0, 0, 0, 0,
                 1, 0, 0, 0, 0),
               nrow = 2, byrow = T)
colnames(wgt_1) = colnames(sampdata_1r)

wgt_1
> wgt_1
     age hgt wgt bmi hc
[1,]   1   0   0   0  0
[2,]   1   0   0   0  0

これらの設定を使ってamputeを実行します。欠測レコードの割合はprop = 0.4としました。

amp_1 = ampute(sampdata_1r, patterns = ptn_1, freq = frq_1, weights = wgt_1, prop = 0.4)
miss_1 = amp_1$amp

意図した通りに欠測データが作成されたかを確認してみましょう。
欠測（NA）の数を見るにはsummaryが手っ取り早いです。

summary(miss_1)

> summary(miss_1)
      age              hgt              wgt              bmi              hc       
 Min.   : 0.035   Min.   : 50.00   Min.   :  3.14   Min.   :11.77   Min.   :33.70  
 1st Qu.: 1.534   1st Qu.: 74.28   1st Qu.: 10.77   1st Qu.:15.86   1st Qu.:48.27  
 Median :10.410   Median :101.40   Median : 30.95   Median :17.37   Median :53.20  
 Mean   : 9.109   Mean   :114.09   Mean   : 35.09   Mean   :17.99   Mean   :51.60  
 3rd Qu.:15.169   3rd Qu.:156.50   3rd Qu.: 56.90   3rd Qu.:19.40   3rd Qu.:56.00  
 Max.   :20.813   Max.   :196.70   Max.   :117.40   Max.   :31.74   Max.   :65.00  
                  NA's   :282      NA's   :92

レコード数684なので、hgtは282 / 684 * 100 = 41.2%、wgtは92 / 684 * 100 = 13.5%が欠測しています。
全体の欠測レコードの割合は40%になるようにしていましたので、hgtは指定に近い欠測割合となっています。
またwgtが欠測する割合は0.4 * 0.3 * 100 = 12%になるよう設定していましたが、こちらもほぼ狙い通りになっているようです。

年齢が高いほど欠測が多いかどうかは、散布図を描いて確認してみます。
元のデータと欠測を発生させたデータのそれぞれについて、横軸を年齢、縦軸を身長とした散布図を作成します。

ggplot(data = miss_1, aes(x = age, y = hgt)) + geom_point() + labs(title = "amputed")
ggplot(data = sampdata_1, aes(x = age, y = hgt)) + geom_point() + labs(title = "original")

f:id:mstour:20210824071114p:plain
f:id:mstour:20210824071124p:plain

上が欠測データ、下が元データです。右の方に行くほどデータ点がスカスカになっているのがわかるでしょうか？

事例2. sleepstudy データセット

sleepstudyデータセットについても、同じ要領で欠測を発生させてみます。
こちらは、毎日繰り返し測定した値が横に並ぶように事前に加工をしていました。
以下のような欠測パターンを考えることにします。

5日目(Day5)から9日目(Day9)までのデータが欠測しうる
ある時点のデータが欠測した場合、以降の時点も欠測する*4
ある時点のデータが欠測するかどうかは、前の時点の値だけに影響される

このデータは変数が多いため、amputeに与えるパラメータの行列を手入力で作るのは少し面倒です。
そこで、デフォルト設定の行列を加工して使用することにします。
まずは何も指定せずにamputeを実行します。

defo_2 = ampute(sampdata_2r)

> defo_2 = ampute(sampdata_2r)
 警告メッセージ: 
Data is made numeric because the calculation of weights requires numeric data

numeric型以外の変数がデータセットに入っていると警告が出ます。
今回は"Subject"（個人の識別番号）がfactor型だったのが原因*5で、amputeされたデータではnumericに変換されていました。
"Subject"は今回の欠測の発生には関与させない予定なので、ここでは無視して進めます。

欠測パターンと重みの設定をそれぞれ抽出し、編集します。
freqについては、直接入力してもいいでしょう。

# 欠測パターン
ptn_defo_2 = defo_2$patterns
ptn_2 = ptn_defo_2[7:11,]
ptn_2[1, 8:11] = c(0, 0, 0, 0)
ptn_2[2, 9:11] = c(0, 0, 0)
ptn_2[3, 10:11] = c(0, 0)
ptn_2[4, 11] = 0

frq_2 = c(0.2, 0.2, 0.2, 0.2, 0.2)

wgt_defo_2 = defo_2$weights
wgt_2 = wgt_defo_2[7:11,]
wgt_2[1, ] = c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)
wgt_2[2, ] = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0)
wgt_2[3, ] = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0)
wgt_2[4, ] = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0)
wgt_2[5, ] = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0)

以下のような行列ができます。

> ptn_2
   Subject Day0 Day1 Day2 Day3 Day4 Day5 Day6 Day7 Day8 Day9
7        1    1    1    1    1    1    0    0    0    0    0
8        1    1    1    1    1    1    1    0    0    0    0
9        1    1    1    1    1    1    1    1    0    0    0
10       1    1    1    1    1    1    1    1    1    0    0
11       1    1    1    1    1    1    1    1    1    1    0

> wgt_2
   Subject Day0 Day1 Day2 Day3 Day4 Day5 Day6 Day7 Day8 Day9
7        0    0    0    0    0    1    0    0    0    0    0
8        0    0    0    0    0    0    1    0    0    0    0
9        0    0    0    0    0    0    0    1    0    0    0
10       0    0    0    0    0    0    0    0    1    0    0
11       0    0    0    0    0    0    0    0    0    1    0

それでは、これらを使ってamputeを実行しましょう。

amp_2 = ampute(sampdata_2r, patterns = ptn_2, freq = frq_2, weights = wgt_2, prop = 0.4)
miss_2 = amp_2$amp

どのようなデータセットが作られたかを確認するために、今度は個人ごとの折れ線グラフを作ってみます。
ggplotで折れ線グラフを作るには、縦型のデータを用意する必要があるので、再度データを加工していきます。

# 縦型に変更
miss_2_l = miss_2 %>%
  pivot_longer(cols = c(-Subject), names_to = "Day", names_prefix = "Day", values_to = "Val") %>%
  mutate(Subject = factor(Subject))

# 折れ線グラフ
ggplot(data = miss_2_l, aes(x = Day, y = Val, group = Subject, color = Subject)) + 
  geom_line() + geom_point() + labs(title = "amputed")

以下のようなグラフができました。
f:id:mstour:20210825201147p:plain