主要アウトカムを24週後のHbA1cとした試験で、最終測定に来なかった患者が一定数出たとすると、来院が途切れた理由が血糖コントロールの悪化そのものにある場合、観測されたデータだけを眺めても、その偏りを確かめる材料はデータの中に残っていません。多重代入で通常置くMARの仮定は、この場面で前提を失います。
そこで、Little（1993）が定式化したパターン混合モデルは、欠測のパターンごとにアウトカムの分布を分けて置き、全体をその混合として書き直します。分けて置くと、データから決まる部分と、研究者が仮定として外から与える部分が、式の上ではっきり切り分けられます。
パターン混合モデルが答える対象
この記事で扱うのは、欠測した患者のアウトカム分布へどのような仮定を置くか、そしてその仮定を動かした時に結論がどこで反転するか、という判断です。MARを仮定した多重代入の手順そのものはMICEによる多重代入：MARの仮定と連鎖方程式で扱っているため、ここでは仮定を外したあとの話に絞ります。
この判断を式で追うために、記号を先に決めておきます。$\mathbf{y}$は1人の対象者について得られる測定値をまとめたベクトルで、24週の測定が1つだけの今回はその値を$Y$と書き、$\mathbf{x}$はベースライン値を含む共変量、$T$は治療群なら1、対照群なら0をとる割付変数、$R$は最終測定が観測された患者で0、欠測した患者で1をとる欠測パターンの指標とします。
パターンごとに分けた同時分布
欠測を含むデータを扱う時は、アウトカムの分布と欠測の起こり方を、ひとまとめの同時分布として書くところから始まります。パターン混合モデルは、この同時分布を欠測のパターンで層別した形へ分解します。
$$f(\mathbf{y}, r \mid \mathbf{x}) = f(\mathbf{y} \mid R = r, \mathbf{x})\,\Pr(R = r \mid \mathbf{x}) \tag{1}$$
式(1)は、アウトカムと欠測パターンの同時分布を、そのパターンに属する患者の中でのアウトカムの条件付き分布と、そのパターンが起きる確率との積へ書き直したものです。$f$は確率密度、$\Pr$は確率を表し、$r$は$R$がとる値です。
この分解が効いてくるのは、右辺の第1項に$R=1$、つまり欠測した患者の分布が、独立した対象として現れる点です。
そのため、データが持っている情報は$R=0$の分布と欠測がどれだけ起きたかという割合だけにとどまり、残りは研究者が置くことになりますが、パターン混合モデルの価値は推定精度ではなく、置いた仮定が式のどこに入ったかが見える点にあります。
選択モデルとの分かれ目
同じ同時分布は、別の順序でも分解できます。欠測の起こりやすさをアウトカムの関数として書く選択モデルは、次の形をとります。
$$f(\mathbf{y}, r \mid \mathbf{x}) = f(\mathbf{y} \mid \mathbf{x})\,\Pr(R = r \mid \mathbf{y}, \mathbf{x}) \tag{2}$$
式(2)は、まず全員に共通のアウトカム分布を置き、そのうえで観測されるかどうかの確率をアウトカムの値に依存させる形で、Heckman型の選択モデルがこの系統に入ります。なお、アウトカムと欠測過程を潜在的なランダム効果で結ぶ共有パラメータモデルは、式(1)とも式(2)とも異なる第三の分解にあたり、この記事では名前を挙げるにとどめます。
両者は同じ同時分布を別の順序で書いたものですが、仮定をどこへ置くかは大きく異なります。選択モデルでは、欠測の確率が未観測のアウトカムへどう依存するかを係数の形で指定するため、その係数の値が臨床的に何を意味するのかを説明しにくくなります。
一方、パターン混合モデルでは、仮定は「中止した患者のHbA1cは、残った患者より0.5%高い水準で推移する」という、アウトカムの単位を持つ量として置かれます。
データから決まらない部分の埋め方
仮定を置く場所がはっきりしたところで、では何を置けばよいのかが次の問題になります。$R=1$の分布は識別できないため、外部から制約を与えないと推定が進みません。
この制約の置き方を、Little（1994）は正規分布の場合について整理しました。最も単純な制約は欠測群の条件付き分布を観測群と同じと置くもので、この制約のもとでの推定値は、共変量を揃えたうえでMARを仮定した多重代入と一致します。
これに対して、臨床的により意味を持つのは、治療を中止した患者が対照群の経過をたどると置く制約です。Carpenterら（2013）は、参照群の分布を中止後へ持ち込むjump to referenceやcopy increments in referenceといった制約を、多重代入の枠組みで実装できる形に整理しました。
この制約も実装は多重代入になるため、統合はMICEによる多重代入：Rでの実行とRubinの統合則と同じ規則で進みます。
臨床的な仮定とズレの量の対応
参照群を使う制約が使えない場面もあります。対照群も同じように悪化している場合や単群試験の場合がそれにあたり、そこで実務でよく使われるのが、観測群から持ち込んだ分布を一定量だけ悪い側へ移す制約です。
$$\mathbb{E}[Y \mid R = 1, T = 1, \mathbf{x}] = \mathbb{E}[Y \mid R = 0, T = 1, \mathbf{x}] + \delta \tag{3}$$
式(3)は、治療群で欠測した患者の平均アウトカムを、同じ治療群にいて同じ共変量を持つ観測された患者の平均より$\delta$だけ高い位置に置くという制約です。$\delta$はアウトカムと同じ単位を持ち、HbA1cを扱う今回の設定では%単位の量になります。
この$\delta$は推定される量ではなく、研究者が値を置く量です。そのため、感度分析としては1つの値を選ぶのではなく、臨床的にありうる範囲を動かして推定値がどう動くかを見ることになります。
具体的には、$\delta=0$から始めて、治療群の中止者が観測された患者よりどれだけ悪い値をとったら結論が変わるかを探す、という読み方が中心になります。
ランダムに生成したデータでの確認
ここまでの制約の置き方が推定値をどちらへ動かすかを、乱数で生成したデータで確かめます。対照群300例、治療群300例の2群並行試験を想定し、ベースラインのHbA1cは平均8.6%、標準偏差0.8%、24週時点の真の治療効果は-0.50%と置きました。
この設定では、脱落は24週のHbA1cが高い患者ほど起きやすく、治療群ではやや起きにくくしてあるため、欠測はMNARにあたります。使用したパッケージはbaseRとstatsだけです。
set.seed(42)
treatment <- rep(c(0, 1), each = 300)
baseline_hba1c <- rnorm(600, mean = 8.6, sd = 0.8)
outcome_true <- 2.20 + 0.65 * baseline_hba1c - 0.50 * treatment +
rnorm(600, mean = 0, sd = 0.60)
dropout_logit <- -1.50 + 0.55 * (baseline_hba1c - 8.6) +
0.80 * (outcome_true - 7.6) - 0.35 * treatment
is_missing <- rbinom(600, size = 1, prob = 1 / (1 + exp(-dropout_logit))) == 1
outcome_observed <- ifelse(is_missing, NA, outcome_true)続いて、欠測の割合と、比較の基準になる2つの治療効果を求めます。全例が観測できた場合の値は実務では手に入らない基準で、これに観測された患者だけを使った完全ケース解析の値を並べます。
fit_full <- lm(outcome_true ~ treatment + baseline_hba1c)
fit_cc <- lm(outcome_observed ~ treatment + baseline_hba1c, subset = !is_missing)
cat(round(c(mean(is_missing), tapply(is_missing, treatment, mean)), 3), "|",
round(c(coef(fit_full)[2], coef(fit_cc)[2]), 4))#> 0.178 0.243 0.113 | -0.4622 -0.4678出力の先頭3つが欠測の割合で、全体では107例、17.8%が欠測し、対照群では24.3%、治療群では11.3%でした。欠測割合が群間で大きく異なるため、観測されたデータだけを使った比較は、両群で異なる大きさの偏りを含みます。
続く2つが治療効果で、全例では-0.4622%、完全ケース解析では-0.4678%でした。差は-0.0056%で、完全ケース解析の標準誤差0.0543%の10分の1ほどにとどまるため、この1回のデータからは偏りの向きも大きさも読み取れません。
読み違えやすいのはここで、2つの値が近かったことは完全ケース解析が安全だったことの確認にはなりません。欠測をMNARに設定しても、1つの乱数系列で得た差は偶然の変動に埋もれるため、偏りの向きを論じるには設定を変えた反復が必要になります。
制約を動かした時の推定値の変化
次に、式(3)の制約を実際に置きます。観測群の回帰から事後分布を引いて欠測値を100組作り、そのうち治療群の中止者にあたる値へ$\delta$を加えてから、各組で同じ回帰を当てはめてRubinの規則で統合しました。
design_all <- cbind(1, treatment, baseline_hba1c)
design_obs <- design_all[!is_missing, ]; y_obs <- outcome_observed[!is_missing]
xtx_inv_obs <- solve(crossprod(design_obs))
beta_hat_obs <- xtx_inv_obs %*% crossprod(design_obs, y_obs)
rss_obs <- sum((y_obs - design_obs %*% beta_hat_obs)^2)
imputed_matrix <- matrix(NA_real_, nrow = sum(is_missing), ncol = 100)
for (m in 1:100) {
sigma2_draw <- rss_obs / rchisq(1, df = nrow(design_obs) - 3)
beta_draw <- beta_hat_obs + t(chol(sigma2_draw * xtx_inv_obs)) %*% rnorm(3)
imputed_matrix[, m] <- design_all[is_missing, ] %*% beta_draw +
rnorm(sum(is_missing), mean = 0, sd = sqrt(sigma2_draw))
}shift_indicator <- as.numeric(is_missing & treatment == 1)
pool_rubin <- function(delta) {
estimates <- numeric(100); variances <- numeric(100)
for (m in 1:100) {
completed <- outcome_observed
completed[is_missing] <- imputed_matrix[, m] + delta * shift_indicator[is_missing]
coefs <- coef(summary(lm(completed ~ treatment + baseline_hba1c)))
estimates[m] <- coefs[2, 1]; variances[m] <- coefs[2, 2]^2
}
total_var <- mean(variances) + 1.01 * var(estimates)
df_rubin <- 99 * (1 + mean(variances) / (1.01 * var(estimates)))^2
mean(estimates) + c(0, -1, 1) * qt(0.975, df_rubin) * sqrt(total_var)
}
cat(round(pool_rubin(0), 4), "|", round(pool_rubin(2.9), 4))#> -0.4675 -0.5745 -0.3604 | -0.1381 -0.286 0.0097出力の前半が$\delta=0$、つまり欠測群の分布を観測群と同じと置いた場合で、推定値は-0.4675%、95%信頼区間は-0.5745から-0.3604でした。この制約はMARを仮定した多重代入と同じ内容のため、パターン混合モデルの枠組みを使った時点でMNARを扱ったことになるわけではありません。
一方、後半の$\delta$を2.9%まで上げた場合は推定値が-0.1381%へ縮み、信頼区間の上端が0.0097となって0をまたぎました。$\delta$を0から3.0%まで0.1%刻みで動かし、同じ統合を繰り返した結果が図1です。
ここで置いた設定について、一点注意が必要です。$\delta$を加えたのは治療群の中止者34例だけで、対照群の中止者73例は観測群と同じ分布のまま残してあり、これは治療に不利な側だけへ振る保守的な感度分析にあたります。
そのため、転換点の値はこの非対称な設定に依存します。対照群の中止者にも同じ量を与えれば両群の平均が同じ向きへ動いて係数の変化はほぼ打ち消され、転換点はまったく異なる値になります。
転換点の臨床的な読み方
図1で読むべきなのは推定値の線そのものより、信頼区間の帯が0と交わるまでの距離です。今回のデータでは、治療群で中止した34例が、同じ治療群で同じベースライン値を持つ観測された患者より24週HbA1cが2.9%高い、という水準まで悪く仮定して、ようやく治療効果の信頼区間が0をまたぎました。
この2.9%という差の大きさを、元のデータの散らばりと比べておきます。ベースラインHbA1cの標準偏差が0.79%でしたから、転換点は標準偏差のおよそ3.7倍にあたります。
ここから先、標準偏差の3.7倍という開きが実際にありうるかどうかは、対象疾患で中止した患者の経過として報告されている値に照らして、研究者が判断する事項です。統計の側から言えるのは、転換点が観測された散らばりに対してどの位置にあるかまでで、その先の判断は領域の知見に委ねられます。
逆に、転換点が測定の散らばりに比べて小さく出る研究では、主解析の結論は欠測した患者の扱い方ひとつで入れ替わるため、点推定値だけを本文に書いて感度分析を補足資料へ回す書き方は適切ではありません。転換点の値を本文へ書き、それが臨床的にありうる差かどうかを考察で評価する方が、読者は結果の頑健さを自分で判断できます。
なお、今回観測された推定値と$\delta$の直線関係は、動かす対象を治療群の中止者に限り、線形回帰を当てはめたことから来ています。別のモデルや別の制約では、この関係は直線になりません。
公的な勧告と報告での扱い
感度分析としての位置づけは、統計家の好みの問題ではなく、公的な勧告にも書かれています。米国National Research Councilの報告を要約したLittleら（2012）は、主解析が置いた欠測の仮定を外した時に結論がどう変わるかを調べる感度分析を、臨床試験の標準的な手続きとして勧告しました。
この勧告に応える具体的な手順を示したのがRatitchら（2013）で、中止後に治療効果が失われる、中止後に悪化する、といった臨床的な想定を、どの制約へ対応させるかを対応表の形で整理しています。パターン混合モデルによる解析を報告する時は、置いた制約の名前、$\delta$を動かした範囲とその範囲の決め方、そして転換点の値まで書くことになります。
なお、MNARを前提とした感度分析の記載事項全般はMICEによる多重代入：MNARの感度分析と論文への記載で扱っているため、ここではパターン混合モデル固有の記載にとどめます。
研究計画で先に決めること
パターン混合モデルを使うかどうかを決める前に、決めておくべきことがあります。欠測した患者について、どのような経過をたどったと考えるのか、その想定を文として書けるかどうかです。
書けるなら、その文はほぼそのまま制約になります。
その想定が書けたうえで、私なら、置いた$\delta$の値が結論を決めてしまいデータが持っている情報との区別がつかなくなるため、主解析はMARを仮定した多重代入で進め、パターン混合モデルは感度分析として添えます。
ただし、中止の理由が有効性の欠如に偏っていることが事前に分かっている場合は、参照群ベースの制約を主解析に据える設計も選択肢になります。アウトカム自体を補完してよいかという判断はMICEによる多重代入：アウトカムを補完してよいかの整理が先に来ます。
また、統計解析計画を書く段階で$\delta$の範囲まで決めておくと、結果を見てから範囲を選んだのではないかという疑いを避けられます。範囲の根拠は、先行研究で報告された中止者のアウトカム、あるいは臨床的に意味のある差の大きさから取ります。
ここまで読んで、自分の研究では$\delta$をどこまで動かせばよいのか、と考えた方もいるかもしれません。実際のデータでは、欠測割合の群間差、中止の理由、主要アウトカムの測定間隔、そして事前に定めた解析計画まで含めて判断することになるため、Dr.データサイエンスでは、こうした欠測データの解析方針についてのご相談を承っています。
参考文献
- Little RJA. Pattern-mixture models for multivariate incomplete data. Journal of the American Statistical Association. 1993;88(421):125-134.
- Little RJA. A class of pattern-mixture models for normal incomplete data. Biometrika. 1994;81(3):471-483.
- Ratitch B, O’Kelly M, Tosiello R. Missing data in clinical trials: from clinical assumptions to statistical analysis using pattern mixture models. Pharmaceutical Statistics. 2013;12(6):337-347.
- Carpenter JR, Roger JH, Kenward MG. Analysis of longitudinal trials with protocol deviation: a framework for relevant, accessible assumptions, and inference via multiple imputation. Journal of Biopharmaceutical Statistics. 2013;23(6):1352-1371.
- Little RJ, D’Agostino R, Cohen ML, et al. The prevention and treatment of missing data in clinical trials. New England Journal of Medicine. 2012;367(14):1355-1360.
