欠測データ分析シミュレーションのための実戦練習:(2)反復測定データの欠測(最終時点のみ)

はじめに

今回は、以前の記事でご紹介したシミュレーションプログラム作成練習の続きをやっていきたいと思います。
mstour.hatenablog.com

題材は、第1回と同様こちらの論文で実施されたシミュレーション研究です。
www.ncbi.nlm.nih.gov

前回は1つのアウトカム変数だけが欠測する状況を扱ったのですが、アウトカム変数が同じ個人から繰り返し測定される状況を考えたいこともあります。少し難易度が上がりますが、現実の場面では役に立つことがあるかと思います。
今回は、アウトカム変数が同じ個人から2回測定され、関心のある最後の測定値が欠測する場合を考えてみます。

データ

今回も、対照治療に対する試験治療の効果を評価するランダム化比較試験を想定します。
まずは完全データの作成手順について、Rコードとともに説明します。

被験者数は2つのグループ(一方は試験治療を、他方は対照治療を受けます)合わせて600人で、ランダム化前に測定された連続型ベースライン変数の1つである X_i*1がアウトカムに影響するという状況を想定します*2。またグループを示す変数を T_iとします。
まず、 X_i T_iを発生させます。 X_iは前回と同様、標準正規分布に従うものとします*3

n = 600
T = c(rep(1, n/2), rep(0, n/2))
X = rnorm(n, mean = 0, sd = 1)

次に、アウトカム変数のうち1回目に測定されたものを Z_i、2回目に測定されたものを Y_iとします。同じ個人から測定される値どうしの相関を考慮するため、 Z_i, Y_iは多変量正規分布からサンプリングします。
また、 Z_i, Y_iの平均はそれぞれ治療グループ T_iごと、ベースライン変数 X_iの値ごとに異なるものとします*4
今回のシミュレーションでは、具体的には以下のように設定します。

 \displaystyle
\left( \begin{array}{c}
Z_i \\ Y_i
\end{array} \right)
\sim
N \left\{ \left( \begin{array}{c}
0.3 T_i + 0.3 X_i \\ 0.3 T_i + 0.3 X_i
\end{array} \right),
\left( \begin{array}{cc}
1 & \sigma_{ZY} \\ \sigma_{ZY} & 1
\end{array} \right)
\right\}

アウトカム変数どうしの相関 \sigma_{ZY}は、0.3と0.7の2パターンの設定を考えることにします。
今回は次のようなコードで作成してみます。
まずは平均0・分散1で、相関が0.3(または0.7)の2変量正規分布からのサンプリングを行います。

# mvrnorm関数を使うため、MASSを読み込みます
library(MASS)

s_z = 1
s_y = 1

# 相関0.3の場合
s_zy_1 = 0.3
cov_zy_1 = matrix(c(s_z, s_zy_1, s_zy_1, s_y), nrow = 2)
e_zy_1 = mvrnorm(n, c(0, 0), cov_zy_1)

# 相関0.7の場合
s_zy_2 = 0.7
cov_zy_2 = matrix(c(s_z, s_zy_2, s_zy_2, s_y), nrow = 2)
e_zy_2 = mvrnorm(n, c(0, 0), cov_zy_2)

続いて、最初に発生させた T_i, X_iを用いて平均を定義し、上記で発生させたサンプルに加えることで Z_i, Y_iを得ます。

alpha1 = 0.3
alpha2 = 0.3
beta1 = 0.3
beta2 = 0.3

# 相関0.3の場合
Z_1 = alpha1 * T + alpha2 * X + e_zy_1[, 1]
Y_1 = beta1 * T + beta2 * X + e_zy_1[, 2]
cd_1 = data.frame(ID = 1:n, T = T, X = X, Z = Z_1, Y = Y_1)

# 相関0.7の場合
Z_2 = alpha1 * T + alpha2 * X + e_zy_2[, 1]
Y_2 = beta1 * T + beta2 * X + e_zy_2[, 2]
cd_2 = data.frame(ID = 1:n, T = T, X = X, Z = Z_2, Y = Y_2)

このような完全データができあがります。

f:id:mstour:20211010194658j:plain

発生させた完全データのうち、今回は Y_iのみが欠測するものとします。
欠測メカニズムはMissing At Random(MAR)とし、欠測確率に関係する変数は前の時点におけるアウトカム変数である Z_iだけという設定にします。 Z_iの影響度は「値が1標準偏差大きくなるごとに欠測のオッズが2.5倍になる」ものとし、欠測するデータの割合は50%とします。

f:id:mstour:20211010175127j:plain

よって、各個人の Y_iの欠測確率は次のように定義できます:

 \displaystyle
\log \left( \frac{\pi_i}{1-\pi_i} \right) = -0.375 + 2.5 Z_i

切片の-0.375は、前回と同様、 Z_iに平均値 0.3 \times 0.5 + 0.3 \times 0 = 0.15を、 \pi_iに0.5を代入して算出した値です。今回もこの設定で欠測データの割合が50%近くなることを確認しました。

欠測を発生させるコード例は以下の通りです。ここではtidyverseパッケージを使っていますので、事前に読み込んでおく必要があります。

library(tidyverse)

b = -0.375

# 相関0.3の場合
missp_1 = 1 / (1 + exp( - (b + 2.5 * Z_1)))
unif_1 = runif(n, min = 0, max = 1)
md_1 = cd_1 %>%
  mutate(missp = missp_1, unif = unif_1) %>%
  mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0), Y_miss = case_when(missfl == 0 ~ Y))

# 相関0.7の場合
missp_2 = 1 / (1 + exp( - (b + 2.5 * Z_2)))
unif_2 = runif(n, min = 0, max = 1)
md_2 = cd_2 %>%
  mutate(missp = missp_2, unif = unif_2) %>%
  mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0), Y_miss = case_when(missfl == 0 ~ Y))

以下のデータの変数「Y_miss」が、 Y_iのうち一部を欠測させたものです。

f:id:mstour:20211010194812j:plain

分析

作成した欠測データに対して、今回は以下の3種類の分析を行います。
推定するのはアウトカム Y_iの2グループ間における平均的な差(治療効果)です。 Y_iの平均は 0.3 T_i + 0.3 X_iなので、推定結果が0.3に近ければ*5適切な分析方法だと考えてよいでしょう。

(1) 欠測を無視した分析(CCA; Complete Case Analysis)
欠測のないレコードのみを使用し、 Y_iを目的変数、 T_i X_iを説明変数とする線形モデルをあてはめます。
以下のようなRコードを用いました。

# モデルあてはめ
mod_2_1 = lm(Y_miss ~ T + X, data = md_2)
# 推定値の抽出
mod_2_1$coefficients["T"]

線形モデルにおける T_iの回帰係数がグループ間のアウトカムの差に相当するので、モデルのあてはめ結果から抽出します。

(2) 多重代入法(MI; Multiple Imputation)
欠測している Y_iを多重代入法により補完した上で、 Y_iのグループ間の単純比較を行います( Y_iを目的変数、 T_iのみを説明変数とする線形モデル)。
 Y_iを補完するモデルの説明変数には、 Y_iの値そのものに影響する T_i X_i、また欠測確率に影響する Z_iを使用します。
Rコードは前回とほぼ同様です。

library(mice)

# 補完モデルを行列形式で設定
pm_2_2 = matrix(c(rep(0, 9), rep(0, 9), rep(0, 9), rep(0, 9),
                  rep(0, 9), rep(0, 9), rep(0, 9), rep(0, 9),
                  0, 1, 1, 1, 0, 0, 0, 0, 0), nrow = 9, byrow = T)
# MICEにより補完
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)[2, "estimate"]

(3) 線形混合モデル(Linear Mixed Model)
反復測定データの分析によく用いられる方法で、欠測がある場合にも(ある条件が満たされていれば)正しい推定を行えることが知られています。(以前書いた説明はこちら
(2)のようなデータの補完はせずに欠測のないレコードのみ使用しますが、(1)と異なり反復測定されたアウトカム変数全て(今回は Y_iおよび Z_i)を目的変数とします。
説明変数は T_i, X_iのほか、測定時点(1、2)を示す情報も含めます。またグループ T_iと測定時点との交互作用項も含めることにします*6

題材の論文では、アウトカム変数の誤差項の共分散行列は「unstructured covariance matrix」とし、自由度の推定にはKenward-Roger法を用いた、との記載があります。
(共分散行列についてはこちらの記事で、Kenward-Roger法についてはこちらの記事でそれぞれ紹介しました。)
しかしながら、Rでは現時点ではunstructured covariance matrixとKenward-Roger法を合わせて使用する方法がないようです(参考文献[2])。
今回のシミュレーションでは自由度計算を使う必要はないので、unstructured covariance matrixの指定のみ行うことにします。参考文献[2]に倣い、nlmeパッケージのgls関数を使用します。
そのために、まずは反復測定データ Z_i, Y_iを横並びの状態から縦型データへ加工します。R以外の統計ソフトウェアでも、混合モデルを用いる際はこの形のデータを必要とすることがよくあります。

# データセットmd_2を縦型に変更
md_2_l = md_2 %>%
  dplyr::select(-Y, -missp, -unif, -missfl) %>%
  tidyr::pivot_longer(cols = c(Z, Y_miss), names_to = "VARS", values_to = "VAL") %>%
  dplyr::mutate(TIME = case_when(VARS == "Z" ~ 1, VARS == "Y_miss" ~ 2))

最初にdplyr::selectで不要な変数を削除し、次にtidyr::pivot_longerで横型のデータセットを縦型に変更しています。cols = で縦に並べる変数を指定します。names_to = で縦型に変更した元の変数名を格納する変数の名称、values_to = で縦型に変更した値を格納する変数の名称をそれぞれ決めることができます。
最後にdplyr::mutateで時点を示す変数を作りました。
縦型に加工した後のデータは次のようになります。

f:id:mstour:20211010194926j:plain

混合モデルはgls関数で次のように実行できます。

library(nlme)

# 誤差項の共分散行列をunstructuredとするLinear mixed modelのあてはめ
mod_2_3 = gls(VAL ~ T + X + factor(TIME) + T*factor(TIME), data = md_2_l, na.action = na.omit,
              correlation = corSymm(form = ~ TIME | ID), weights = varIdent(form = ~ 1 | TIME))
# 推定結果の抽出(時点2の結果を見るため、交互作用項の係数を足す)
mod_2_3$coefficients["T"] + mod_2_3$coefficients["T:factor(TIME)2"]

上記のコードでは、correlation = の部分で共分散行列の非対角成分、weights = の部分で対角成分をそれぞれ設定しています。またna.action = na.omit と指定することで、欠測のある個人が分析対象から除外されてしまわないようにしています*7

今回のシミュレーションでも、以上の一連の流れを所定の回数反復し、推定値の平均と標準偏差をプロットすることにしました。
コード全体は次項にまとめています(作図のコードもこちらに掲載しています)。

シミュレーションコードまとめ

library(tidyverse)
library(mice)
library(MASS)
library(nlme)

set.seed(1815)

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

n = 600
alpha1 = 0.3
alpha2 = 0.3
beta1 = 0.3
beta2 = 0.3
s_z = 1
s_y = 1
s_zy_1 = 0.3
s_zy_2 = 0.7
cov_zy_1 = matrix(c(s_z, s_zy_1, s_zy_1, s_y), nrow = 2)
cov_zy_2 = matrix(c(s_z, s_zy_2, s_zy_2, s_y), nrow = 2)
b = -0.375

### corr=0.3
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_zy_1 = mvrnorm(n, c(0, 0), cov_zy_1)
  Z_1 = alpha1 * T + alpha2 * X + e_zy_1[, 1]
  Y_1 = beta1 * T + beta2 * X + e_zy_1[, 2]
  cd_1 = data.frame(ID = 1:n, T = T, X = X, Z = Z_1, Y = Y_1)
  missp_1 = 1 / (1 + exp( - (b + 2.5 * Z_1)))
  unif_1 = runif(n, min = 0, max = 1)
  md_1 = cd_1 %>%
    mutate(missp = missp_1, unif = unif_1) %>%
    mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0),
           Y_miss = case_when(missfl == 0 ~ Y))
  ### analysis
  # 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, 9), rep(0, 9), rep(0, 9), rep(0, 9),
                    rep(0, 9), rep(0, 9), rep(0, 9), rep(0, 9),
                    0, 1, 1, 1, 0, 0, 0, 0, 0), nrow = 9, 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"])
  # linear mixed model
  md_1_l = md_1 %>%
    dplyr::select(-Y, -missp, -unif, -missfl) %>%
    tidyr::pivot_longer(cols = c(Z, Y_miss), names_to = "VARS", values_to = "VAL") %>%
    dplyr::mutate(TIME = case_when(VARS == "Z" ~ 1, VARS == "Y_miss" ~ 2))
  mod_1_3 = gls(VAL ~ T + X + factor(TIME) + T*factor(TIME),
                data = md_1_l, na.action = na.omit,
                correlation = corSymm(form = ~ TIME | ID),
                weights = varIdent(form = ~ 1 | TIME))
  tmp = mod_1_3$coefficients["T"] + mod_1_3$coefficients["T:factor(TIME)2"]
  result_1_3 = rbind(result_1_3, tmp)
}

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", "Linear mixed model"),
                       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.2, 0.3, 0.4, 0.5), limits = c(0.10, 0.51)) +
  labs(title = "Corr(Y, Z) = 0.30", x = "", y = "Treatment effect estimate") +
  theme_classic()

### corr=0.7
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_zy_2 = mvrnorm(n, c(0, 0), cov_zy_2)
  Z_2 = alpha1 * T + alpha2 * X + e_zy_2[, 1]
  Y_2 = beta1 * T + beta2 * X + e_zy_2[, 2]
  cd_2 = data.frame(ID = 1:n, T = T, X = X, Z = Z_2, Y = Y_2)
  missp_2 = 1 / (1 + exp( - (b + 2.5 * Z_2)))
  unif_2 = runif(n, min = 0, max = 1)
  md_2 = cd_2 %>%
    mutate(missp = missp_2, unif = unif_2) %>%
    mutate(missfl = case_when(missp > unif ~ 1, TRUE ~ 0),
           Y_miss = case_when(missfl == 0 ~ Y))
  ### analysis
  # 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, 9), rep(0, 9), rep(0, 9), rep(0, 9),
                    rep(0, 9), rep(0, 9), rep(0, 9), rep(0, 9),
                    0, 1, 1, 1, 0, 0, 0, 0, 0), nrow = 9, 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"])
  # linear mixed model
  md_2_l = md_2 %>%
    dplyr::select(-Y, -missp, -unif, -missfl) %>%
    tidyr::pivot_longer(cols = c(Z, Y_miss), names_to = "VARS", values_to = "VAL") %>%
    dplyr::mutate(TIME = case_when(VARS == "Z" ~ 1, VARS == "Y_miss" ~ 2))
  mod_2_3 = gls(VAL ~ T + X + factor(TIME) + T*factor(TIME),
                data = md_2_l, na.action = na.omit,
                correlation = corSymm(form = ~ TIME | ID),
                weights = varIdent(form = ~ 1 | TIME))
  tmp = mod_2_3$coefficients["T"] + mod_2_3$coefficients["T:factor(TIME)2"]
  result_2_3 = rbind(result_2_3, tmp)
}

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", "Linear mixed model"),
                       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.2, 0.3, 0.4, 0.5), limits = c(0.10, 0.51)) +
  labs(title = "Corr(Y, Z) = 0.70", x = "", y = "Treatment effect estimate") +
  theme_classic()

結果

今回も、治療効果(グループ間のアウトカムの差)の推定が平均的にどのようになるか、3つの分析方法それぞれの結果をプロットして確認します*8

まず、アウトカム変数の時点間の相関を0.3とした場合です。

f:id:mstour:20211010071135p:plain

題材の論文の結果と同様に、多重代入法(青色)と混合モデル(緑色)はいずれも平均的に真の値0.3をほぼ正しく推定していますが、欠測を無視した分析(赤色)では少し小さめに評価してしまう傾向にありました。

相関が0.7の場合は以下のようになりました。

f:id:mstour:20211010072242p:plain

全シミュレーションを通した平均の値(丸印)は論文と少し異なる結果になりましたが、相関が0.3の場合と比べて欠測を無視したときの過小推定の程度が大きくなっている様子が再現できています。

おわりに

今回は同じ個人からアウトカム変数が2回測定されている反復測定データを考え、最終時点のみが欠測する場合の分析結果のシミュレーションを行いました。
現実には、もっと多くの測定回数があったり、複数時点で欠測が発生するような状況も多いかと思います。
それらも今回と同様の方法でシミュレーションすることができますが、せっかくなので次回以降に実施してみます。題材の論文ではこのようなケースは取り扱われていなかったので、次はお手本なしでいきます。

参考文献

[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] Mixed model repeated measures (MMRM) in Stata, SAS and R – The Stats Geek

*1:添字 iは個人を示す記号です。よって、サンプルサイズがNの場合 X_1から X_Nまでの独立なN個のデータがあると考えてください。その他の変数も同様です。

*2:例えば年齢やBMIといった背景因子や、アウトカム変数そのものの介入前値などが考えられます。

*3:発生させたデータに「 sを掛けて」「 mを足す」という操作をすれば、平均 m、分散 s^2正規分布からのサンプリングと同等になります。

*4:以下の平均の定義で、 T_i, X_iがアウトカムの平均に与える影響をそれぞれ0.3と設定しています。

*5:理論的にはどういう意味で0.3に近いのか(例えば、無限に同じ試験を繰り返した時の平均が0.3に近いのか、サンプルサイズを無限に増やした時に0.3に近くなるのか)を考える必要がありますが、筆者の手に余る話なのでここでは深く考えないことにします。

*6:交互作用を含めることで、時点1ではグループ間に差はないけれど時点2では差が生じている、という状況を表現することができます。

*7:長くなるので、詳細は参考文献[2]をご参照ください。

*8:題材の論文とは横軸の順番と色付けが異なります。うまく調整する方法は調べればあると思うのですが、本題と外れるので今回はそこまで踏み込んでいません。