欠測データ分析シミュレーションのための実戦練習：(1)1変数の欠測

はじめに

過去2回、Rのmiceパッケージを用いた欠測データの作成方法を紹介してきました。
（過去の記事は以下です。）
mstour.hatenablog.com
mstour.hatenablog.com

今回はもう少し踏み込んで、欠測データ分析に関する実際のシミュレーション研究の再現をしてみたいと思います。題材はSullivanらによる下記の論文です。
（こちらはどの環境からでもアクセスできると思うのですが、見れなかったらすみません。）
www.ncbi.nlm.nih.gov

この論文は欠測メカニズムがMARである様々な欠測データを作成し、ランダム化比較試験における多重代入法・混合効果モデル・欠測を無視した解析の性能比較などを行なったものです。

なお、今回は色々事情があってmiceパッケージのampute関数は使わないことにしました。
理由は、各変数の欠測への影響度を設定することが思いのほか難しいためです。

どういうことかと言うと、複数の変数間の相対的な重みの設定は引数weightsにより行えるのですが、ある変数の絶対的な影響度、つまり回帰係数を具体的にいくつにするという設定を直接行う方法がないようなのです。具体的には変数 $y_i$ の欠測確率を以下のロジスティック回帰モデル

$\displaystyle \log \left( \frac{\pi_i}{1-\pi_i} \right) = \beta_0 + \beta_1 x_i + \beta_2 z_i$

で与えるものとすると、変数 $x_i$ と $z_i$ のどちらがより欠測への影響が大きいかは引数weightsで指定することができる一方、係数 $\beta_1$ と $\beta_2$ の値を直接指定することができないということです。

もちろん上手く設定する方法はあるのかもしれませんが・・・色々な資料やソースコードを見る限りでは現時点で解決することができませんでしたので、今回は自分でプログラムを作って欠測データを発生させることにします。

シミュレーション研究の設定

さて、今回題材とする研究は、ランダム化比較試験でデータの欠測が発生した場合に、多重代入法はお勧めできる対処法なのか？を検討したものです。
本記事では、最初の例として挙げられている以下のシミュレーションを行い、論文内のFigure 1を再現してみます。

2つの治療をランダムに割り当てる状況を想定し、アウトカムにはある一つの共変量が影響しているものとします。
まずは欠測のない完全データを用意します。
具体的にはアウトカムを $Y_i$ 、治療を表す変数を $T_i$ （試験治療の場合に1、対照治療の場合に0となるものとします）、共変量を $X_i$ 、誤差を $e_i$ とし、次のモデルからデータを生成します。

$\displaystyle Y_i = 0.30 T_i + 0.70 X_i + e_i \hspace{20pt} (i = 1, \cdots. N)$

ただし $X_i, e_i$ はそれぞれ互いに独立に標準正規分布に従うものとします。
Rを用いたデータの生成方法は後述します。

次に、完全データのアウトカム $Y_i$ の一部を人為的に欠測させます。
欠測のパターンとして以下の3つを考え、それぞれに対応する欠測データを作ります。

欠測パターン1（MAR X）：欠測メカニズムはMARで、欠測はXにのみ依存する。Xが1標準偏差大きくなるごとに、欠測するオッズ*1は2.5大きくなる。

欠測パターン2（MAR X+T）：欠測メカニズムはMARで、欠測はXおよびTに依存する。Xが1標準偏差大きくなるごとに、欠測するオッズは2.5大きくなる。また対照治療の欠測するオッズは試験治療より2.5大きい。

欠測パターン3（MAR X×T）：欠測メカニズムはMARで、欠測はXおよびTに依存する。Xが1標準偏差大きくなるごとに、欠測するオッズは2.5大きくなる。また $X >0$ のとき対照治療の欠測オッズが試験治療より2.5大きくなり、 $X \leq 0$ のとき試験治療の欠測オッズが対照治療より2.5大きくなる。

ちょっとややこしくなりましたが、ロジスティック回帰モデルで上記1〜3の欠測確率を表現すると、それぞれ次のようになります。

欠測パターン1：

$\displaystyle \log \left( \frac{\pi_i}{1-\pi_i} \right) = \beta_0 + 2.5 X_i \\$

欠測パターン2：

$\displaystyle \log \left( \frac{\pi_i}{1-\pi_i} \right) = \beta_0 + 2.5 X_i - 2.5 T_i \\$

欠測パターン3*2：

$\displaystyle \log \left( \frac{\pi_i}{1-\pi_i} \right) = \beta_0 + 2.5 X_i - 2.5 T_i I(X_i > 0) + 2.5 T_i I(X_i \leq 0)$

$\beta_0$ は、 $Y_i$ が欠測する割合を何%にするかに応じて適切な値を決める必要があります。

今回は欠測の割合を50%に設定するのですが、対応する $\beta_0$ の値を数学的にまともに算出するのはかなり大変なようですので*3、今回はロジスティック回帰モデルの式を $\beta_0$ について整理した上で、 $X_i, T_i$ にそれぞれの期待値を、 $\pi_i$ に0.5を代入するというざっくりした方法で $\beta_0$ を求めます。
具体的には以下の数値になりました。後ほど、欠測割合が意図する通りになっていることを確認します。

欠測パターン1： $\beta_0 = 0$
欠測パターン2： $\beta_0 = 1.25$
欠測パターン3： $\beta_0 = 0$

最後に、作成した欠測データに対して分析を行います。分析方法は以下の3つを試します。

TとXを説明変数とする線形モデル（欠測に関与するXを考慮した分析）
TとXを用いてYを多重代入法で補完し、補完されたデータにTだけを説明変数とする線形モデル*4を使用（欠測に関与するXを考慮した分析）
Tだけを説明変数とする線形モデル（欠測に関与するXを無視した分析）

以上の「完全データ生成」「欠測の発生」「分析」の流れを所定の回数反復することによって、それぞれの方法で平均的に正しい推定が行えるかを評価します。

かなり前置きが長くなりました。では以下で実際にRを使ってシミュレーションをしていきましょう。

シミュレーションコード

シミュレーションを行う前に、まずは毎回行う処理のコードを説明します。
最初に、前述の設定に沿って完全データを作成します。

library(tidyverse)
set.seed(1815)

n = 600
beta1 = 0.3
beta2 = 0.7
T = c(rep(1, n/2), rep(0, n/2))
X = rnorm(n, mean = 0, sd = 1)
e = rnorm(n, mean = 0, sd = 1)
Y = beta1 * T + beta2 * X + e
cd = data.frame(ID = 1:n, T = T, X = X, Y = Y)

試験治療（T=1）と対照治療（T=0）をランダムに割り付けているため、共変量Xの分布が両グループで等しくなるようにデータを生成する必要があります。上記のように完全にランダムに抽出すれば問題ないでしょう。

次に欠測値を発生させます。欠測パターン1〜3に対応する欠測データはそれぞれ以下のように作成できます。

### MAR X
b_1 = 0
missp_1 = 1 / (1 + exp( - (b_1 + 2.5 * X)))
unif_1 = runif(n, min = 0, max = 1)
md_1 = cd %>%
  mutate(missp = missp_1, unif = unif_1) %>%
  mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0),
         Y_miss = case_when(missfl == 0 ~ Y))

### MAR X+T
b_2 = 1.25
missp_2 = 1 / (1 + exp( - (b_2 + 2.5 * X - 2.5 * T)))
unif_2 = runif(n, min = 0, max = 1)
md_2 = cd %>%
  mutate(missp = missp_2, unif = unif_2) %>%
  mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0),
         Y_miss = case_when(missfl == 0 ~ Y))

### MAR X*T
b_3 = 0
XI = ifelse(X>0, 1, 0)
missp_3 = 1 / (1 + exp( - (b_3 + 2.5 * X - 2.5 * T * XI + 2.5 * T * (1-XI))))
unif_3 = runif(n, min = 0, max = 1)
md_3 = cd %>%
  mutate(missp = missp_3, unif = unif_3) %>%
  mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0),
         Y_miss = case_when(missfl == 0 ~ Y))

最初の「MAR X」では、まずmissp_1でロジスティック回帰モデルから各個人の欠測確率を計算し、次に0〜1の値をとる一様乱数を人数分発生させ、乱数の値が欠測確率より小さい場合*5にその個人のデータを欠測させています。他のパターンでも同様です。

欠測割合を50%としていましたが、意図した通りに欠測データが作成できているかを確認しておきましょう。

> sum(is.na(md_1$Y_miss)) / n
[1] 0.4933333

> sum(is.na(md_2$Y_miss)) / n
[1] 0.49

> sum(is.na(md_3$Y_miss)) / n
[1] 0.51

どれも問題ないですね。

欠測データに対する分析は、例えば次のようなコードで実施してみましょう。
ここでは欠測パターン2（XとTが欠測に関与）を例としていますが、他のパターンでも全く同じ方法でOKです。

#----- analysis example(MAR X+T) -----
library(mice)

# adjusted cca
mod_2_1 = lm(Y_miss ~ T + X, data = md_2)
summary(mod_2_1)

# mi overall
pm_2_2 = matrix(c(rep(0, 8), rep(0, 8), rep(0, 8), rep(0, 8),
                  rep(0, 8), rep(0, 8), rep(0, 8),
                  0, 1, 1, 0, 0, 0, 0, 0), nrow = 8, byrow = T)
imp_2_2 = mice(data = md_2, m = 50, method = "norm", predictorMatrix = pm_2_2)
mod_2_2_0 = with(imp_2_2, lm(Y_miss ~ T))
mod_2_2 = pool(mod_2_2_0)
summary(mod_2_2)

# unadjusted cca
mod_2_3 = lm(Y_miss ~ T, data = md_2)
summary(mod_2_3)

多重代入法（上記コードの"mi overall"）は、あらかじめmiceパッケージを読み込んでMICE法で行っています。
欠測の補完モデルの説明変数としてXとTを使う設定にしています。
今回の本題から外れるので詳しい説明は省略し、また別の記事でこの辺りの話を書きたいと思います。
"mod_2_1"、"mod_2_2"、"mod_2_3"に3つの方法それぞれの分析結果を格納しています。

今回実施したシミュレーション用コードの全体は以下の通りです。
所定の繰り返しで得られた治療効果推定値（変数Tの回帰係数です）の平均と標準偏差 *6を計算し、題材の論文と同じようなプロットを作成しました。

#----- simulation -----
### parameter setting
nsim = 2000

n = 600
beta1 = 0.3
beta2 = 0.7

b_1 = 0
b_2 = 1.25
b_3 = 0

### MAR X
result_1_1 = NULL
result_1_2 = NULL
result_1_3 = NULL
for (i in 1:nsim){
  # data
  T = c(rep(1, n/2), rep(0, n/2))
  X = rnorm(n, mean = 0, sd = 1)
  e = rnorm(n, mean = 0, sd = 1)
  Y = beta1 * T + beta2 * X + e
  cd = data.frame(ID = 1:n, T = T, X = X, Y = Y)
  missp_1 = 1 / (1 + exp( - (b_1 + 2.5 * X)))
  unif_1 = runif(n, min = 0, max = 1)
  md_1 = cd %>%
    mutate(missp = missp_1, unif = unif_1) %>%
    mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0),
           Y_miss = case_when(missfl == 0 ~ Y))
  ### analysis
  # adjusted cca
  mod_1_1 = lm(Y_miss ~ T + X, data = md_1)
  result_1_1 = rbind(result_1_1, mod_1_1$coefficients["T"])
  # mi overall
  pm_1_2 = matrix(c(rep(0, 8), rep(0, 8), rep(0, 8), rep(0, 8),
                    rep(0, 8), rep(0, 8), rep(0, 8),
                    0, 1, 1, 0, 0, 0, 0, 0), nrow = 8, byrow = T)
  imp_1_2 = mice(data = md_1, m = 50, method = "norm", predictorMatrix = pm_1_2)
  mod_1_2_0 = with(imp_1_2, lm(Y_miss ~ T))
  mod_1_2 = pool(mod_1_2_0)
  result_1_2 = rbind(result_1_2, summary(mod_1_2)[2, "estimate"])
  # unadjusted cca
  mod_1_3 = lm(Y_miss ~ T, data = md_1)
  result_1_3 = rbind(result_1_3, mod_1_3$coefficients["T"])
}
mean_1 = c(mean(result_1_1), mean(result_1_2), mean(result_1_3))
sd_1 = c(sd(result_1_1), sd(result_1_2), sd(result_1_3))
summary_1 = data.frame(analysis = c("CCA", "MI overall m=50", "Unadjusted CCA"),
                       mean = mean_1, sd = sd_1)
ggplot(data = summary_1, aes(x = analysis, y = mean, color = analysis)) +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd), width = 0.1) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0.3, linetype = 2, size = 0.3) +
  scale_y_continuous(breaks = c(0.0, 0.3, 0.6, 0.9), limits = c(-0.01, 0.91)) +
  labs(title = "MAR X", x = "", y = "Treatment effect estimate") +
  theme_classic()

### MAR X+T
result_2_1 = NULL
result_2_2 = NULL
result_2_3 = NULL
for (i in 1:nsim){
  # data
  T = c(rep(1, n/2), rep(0, n/2))
  X = rnorm(n, mean = 0, sd = 1)
  e = rnorm(n, mean = 0, sd = 1)
  Y = beta1 * T + beta2 * X + e
  cd = data.frame(ID = 1:n, T = T, X = X, Y = Y)
  missp_2 = 1 / (1 + exp( - (b_2 + 2.5 * X - 2.5 * T)))
  unif_2 = runif(n, min = 0, max = 1)
  md_2 = cd %>%
    mutate(missp = missp_2, unif = unif_2) %>%
    mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0),
           Y_miss = case_when(missfl == 0 ~ Y))
  ### analysis
  # adjusted cca
  mod_2_1 = lm(Y_miss ~ T + X, data = md_2)
  result_2_1 = rbind(result_2_1, mod_2_1$coefficients["T"])
  # mi overall
  pm_2_2 = matrix(c(rep(0, 8), rep(0, 8), rep(0, 8), rep(0, 8),
                    rep(0, 8), rep(0, 8), rep(0, 8),
                    0, 1, 1, 0, 0, 0, 0, 0), nrow = 8, byrow = T)
  imp_2_2 = mice(data = md_2, m = 50, method = "norm", predictorMatrix = pm_2_2)
  mod_2_2_0 = with(imp_2_2, lm(Y_miss ~ T))
  mod_2_2 = pool(mod_2_2_0)
  result_2_2 = rbind(result_2_2, summary(mod_2_2)[2, "estimate"])
  # unadjusted cca
  mod_2_3 = lm(Y_miss ~ T, data = md_2)
  result_2_3 = rbind(result_2_3, mod_2_3$coefficients["T"])
}
mean_2 = c(mean(result_2_1), mean(result_2_2), mean(result_2_3))
sd_2 = c(sd(result_2_1), sd(result_2_2), sd(result_2_3))
summary_2 = data.frame(analysis = c("CCA", "MI overall m=50", "Unadjusted CCA"),
                       mean = mean_2, sd = sd_2)
ggplot(data = summary_2, aes(x = analysis, y = mean, color = analysis)) +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd), width = 0.1) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0.3, linetype = 2, size = 0.3) +
  scale_y_continuous(breaks = c(0.0, 0.3, 0.6, 0.9), limits = c(-0.01, 0.91)) +
  labs(title = "MAR X+T", x = "", y = "Treatment effect estimate") +
  theme_classic()

### MAR X*T
result_3_1 = NULL
result_3_2 = NULL
result_3_3 = NULL
for (i in 1:nsim){
  # data
  T = c(rep(1, n/2), rep(0, n/2))
  X = rnorm(n, mean = 0, sd = 1)
  e = rnorm(n, mean = 0, sd = 1)
  Y = beta1 * T + beta2 * X + e
  cd = data.frame(ID = 1:n, T = T, X = X, Y = Y)
  XI = ifelse(X>0, 1, 0)
  missp_3 = 1 / (1 + exp( - (b_3 + 2.5 * X - 2.5 * T * XI + 2.5 * T * (1-XI))))
  unif_3 = runif(n, min = 0, max = 1)
  md_3 = cd %>%
    mutate(missp = missp_3, unif = unif_3) %>%
    mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0),
           Y_miss = case_when(missfl == 0 ~ Y))
  ### analysis
  # adjusted cca
  mod_3_1 = lm(Y_miss ~ T + X, data = md_3)
  result_3_1 = rbind(result_3_1, mod_3_1$coefficients["T"])
  # mi overall
  pm_3_2 = matrix(c(rep(0, 8), rep(0, 8), rep(0, 8), rep(0, 8),
                    rep(0, 8), rep(0, 8), rep(0, 8),
                    0, 1, 1, 0, 0, 0, 0, 0), nrow = 8, byrow = T)
  imp_3_2 = mice(data = md_3, m = 50, method = "norm", predictorMatrix = pm_3_2)
  mod_3_2_0 = with(imp_3_2, lm(Y_miss ~ T))
  mod_3_2 = pool(mod_3_2_0)
  result_3_2 = rbind(result_3_2, summary(mod_3_2)[2, "estimate"])
  # unadjusted cca
  mod_3_3 = lm(Y_miss ~ T, data = md_3)
  result_3_3 = rbind(result_3_3, mod_3_3$coefficients["T"])
}
mean_3 = c(mean(result_3_1), mean(result_3_2), mean(result_3_3))
sd_3 = c(sd(result_3_1), sd(result_3_2), sd(result_3_3))
summary_3 = data.frame(analysis = c("CCA", "MI overall m=50", "Unadjusted CCA"),
                       mean = mean_3, sd = sd_3)
ggplot(data = summary_3, aes(x = analysis, y = mean, color = analysis)) +
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd), width = 0.1) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0.3, linetype = 2, size = 0.3) +
  scale_y_continuous(breaks = c(0.0, 0.3, 0.6, 0.9), limits = c(-0.01, 0.91)) +
  labs(title = "MAR X×T", x = "", y = "Treatment effect estimate") +
  theme_classic()

分析結果

それぞれの欠測パターンについて、以下のような結果になりました。
f:id:mstour:20210917220928p:plain
f:id:mstour:20210917220937p:plain
f:id:mstour:20210917220945p:plain
論文のFigure 1と同じような図になっていますよね？
論文の結果と同じく、欠測に関与する変数を考慮しない"Unadjusted CCA"では推定にバイアスが入ることがわかります*7。
ただしMAR X+Tの場合のみ、論文よりも推定値がかなり大きくなってしまいました（論文では丸印で示す平均が0.45くらいだが、今回の結果は0.6くらい）。
シミュレーション中に偶然大きな値ばかりが出たのかもしれませんし、論文で書かれていない何らかの情報があるのかもしれませんが、ちょっと原因がわかりませんでした。何かわかったら追記したいと思います。

この結果で重要なのは、どの欠測パターンでも、"CCA"（欠測に関与する変数XとT両方を説明変数とする）と"MI overall"（XとTを説明変数とする補完モデルを用いて多重代入法で欠測を補完する）のどちらもバイアスなく変数Tの効果（治療効果）が推定できているという点です。
欠測メカニズムがMARであれば、欠測に関与する変数で条件付けした分析*8を行えばバイアスなく推定できるため、今回の設定においてはわざわざ多重代入法によって欠測値の補完をしなくてもよい、ということになります。
むしろ多重代入法は通常の線形モデルよりも推定のばらつきが大きくなりがちなので、この場合にはむしろ好ましくないと言えるでしょう*9。
今回作成した図でも、"MI overall"の標準偏差の範囲が"CCA"よりも少しだけ広くなっているのがわかるかと思います。

おわりに

今回は欠測データ分析に関するシミュレーションの実践練習として、統計の研究論文の結果を再現することを試みました。
実務や研究では、経時的データの欠測や多変量の欠測など、もっと複雑な状況のシミュレーションを行わないといけないこともあるかと思います。
次回はもう少し難しい問題を取り上げてみます。

参考文献

[1] Sullivan, T. R., White, I. R., Salter, A. B., Ryan, P., & Lee, K. J. (2018). Should multiple imputation be the method of choice for handling missing data in randomized trials?. Statistical methods in medical research, 27(9), 2610-2626.
[2] Van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate imputation by chained equations in R. Journal of statistical software, 45(1), 1-67.
[3] Schouten, R. M., Lugtig, P., & Vink, G. (2018). Generating missing values for simulation purposes: a multivariate amputation procedure. Journal of Statistical Computation and Simulation, 88(15), 2909-2930.

*1:少し分かりづらいですが、ロジスティック回帰モデルで欠測確率を計算するためこのような表現になっています。厳密には異なりますが、ここでは「欠測するオッズは2.5大きくなる」は「欠測確率が2.5倍になる」とだいたい同じようなものだと考えてもらっても構いません。

*2: $I( )$ は、かっこの中身の条件に該当する場合に1を、そうでない場合に0となる関数です。例えば $I(X_i > 0)$ は $X_i$ が正の値のときに1、そうでないとき0をとります。