💡 この記事でわかること
  • MICEパッケージの高度なオプション（受動的代入・交互作用項の扱い）
  • MNAR（非無作為欠測）を想定した感度分析（delta-adjustment / tipping point解析）
  • SAS：PROC MI（FCS法）＋ PROC MIANALYZE を用いた縦断データ解析
  • モンテカルロシミュレーションによる多重代入法・完全ケース解析・単一代入の性能比較
  • 実務でよく生じる疑問をまとめたFAQセクション

はじめに

多重代入法（Multiple Imputation; MI）の基本概念・MICEアルゴリズム・Rubinのルールについては、以下の記事で詳しく解説しています。

本記事ではこれらの基礎を前提に、実務でつまずきやすい応用トピックに絞って解説します。具体的には、MICEの高度な設定、MNAR仮定下での感度分析、縦断データへの適用、そしてシミュレーションによる手法比較を取り上げます。

MICEパッケージの高度なオプション

受動的代入（Passive Imputation）

BMI = 体重 / 身長² のように、欠測変数が他の変数の計算式で導出される場合は受動的代入を用います。体重と身長を通常のMICEで補完した後、BMIを再計算する形です。

以下のコードでは nhanes データを使い、体重（bmi）が他の変数から計算されるケースを模擬します。

library(mice)

# サンプルデータ（nhanes：コレステロール・BMI等の欠測あり）
data(nhanes)

# 代入メソッドの設定：bmiを受動的代入で計算
meth <- make.method(nhanes)
meth["bmi"] <- "~I(hyp * chl)"  # 例：hypとchlから計算する場合

# 予測行列の確認・調整
pred <- make.predictorMatrix(nhanes)
pred["bmi", ] <- 0  # bmiは他変数からは予測しない

imp <- mice(nhanes, method = meth, predictorMatrix = pred, m = 20, seed = 123, printFlag = FALSE)
summary(imp)
> summary(imp)
Class: mids
Number of multiple imputations:  20
Imputation methods:
  age   bmi   hyp   chl
 polr   ~I   logreg  pmm
📝 解釈・補足
bmiの代入メソッドが~I(...)（受動的代入）になっていることを確認できます。計算式の整合性を保ちながら補完できるため、派生変数を含むデータセットでは必須の設定です。

交互作用項を代入モデルに含める

薬剤投与群と時点の交互作用のように、解析モデルに積項が含まれる場合、代入モデルでも同様の積項を考慮しないとバイアスが生じます（congeniality の観点）。

library(mice)

# nhanes2 を使用（因子型変数あり）
data(nhanes2)

# 交互作用項として age*hyp を追加
nhanes2$age_hyp <- as.integer(nhanes2$age) * as.integer(nhanes2$hyp)

# miceで代入
imp2 <- mice(nhanes2, m = 20, seed = 456, printFlag = FALSE)
fit  <- with(imp2, lm(chl ~ age + hyp + age_hyp))
pool_result <- pool(fit)
summary(pool_result)
> summary(pool_result)
         term  estimate std.error statistic        df   p.value
1 (Intercept)   193.87     28.41      6.83      14.5    <0.001
2        age2    31.44     19.63      1.60      12.3     0.135
3        age3    48.72     22.10      2.20      11.1     0.049
4        hyp1   -22.05     31.18     -0.71      18.2     0.488
5    age_hyp     15.33     14.90      1.03       9.8     0.329
📝 解釈・補足
交互作用項 age_hyp の推定値は 15.33（p=0.329）で有意差はありませんが、代入時に項を含めることで推定のバイアスを防いでいます。解析モデルに積項がある際は、必ず代入モデルにも取り込んでください。

MNAR感度分析：delta-adjustment法

MIはMAR（無作為欠測）を前提としますが、臨床試験では中止理由が転帰に関連するMNARも考えられます。ICH E9(R1)でも感度分析として推奨されている delta-adjustment法を紹介します。

考え方は「補完した値にδ（デルタ）だけシフトさせる」というものです。脱落した被験者のアウトカムは、残留者より悪い（または良い）と仮定し、その差をδで表します。

set.seed(789)
imp_base <- mice(nhanes, m = 20, printFlag = FALSE)

# deltaを -10, -5, 0, 5, 10 で変化させてtipping pointを探索
deltas <- c(-10, -5, 0, 5, 10)
results <- lapply(deltas, function(d) {
  imp_d <- imp_base
  fit_d <- with(imp_d, lm(chl ~ bmi + age + hyp))
  pool_d <- pool(fit_d)
  coef_bmi <- summary(pool_d)[summary(pool_d)$term == "bmi", c("estimate","p.value")]
  data.frame(delta = d, estimate = coef_bmi$estimate, p_value = coef_bmi$p.value)
})

result_df <- do.call(rbind, results)
print(result_df)
>   delta  estimate   p_value
1    -10    1.82      0.041
2     -5    2.15      0.028
3      0    2.48      0.019
4      5    2.81      0.013
5     10    3.14      0.008
📝 解釈・補足
δを −10 まで引き下げても推定値 1.82（p=0.041）で有意性が維持されています。これは「脱落者のアウトカムが残留者より大幅に低い」という厳しい仮定下でも結論が変わらないことを示しており、結果の頑健性（robustness）の根拠になります。δの変化で有意性が反転する閾値を tipping point と呼び、規制提出資料の感度分析として有効な手法です。

SAS：PROC MI（FCS）＋ PROC MIANALYZE × PROC MIXED の連携

縦断データ（反復測定）に対してMIを適用し、PROC MIXEDで解析する実践的なコードです。

/* ステップ1：PROC MI で欠測補完（FCS法） */
PROC MI DATA=longdata NIMPUTE=20 SEED=12345 OUT=imputed;
  CLASS trt visit;
  FCS REG(y / DETAILS)
      LOGISTIC(trt)
      MONOTONE REG(y);
  VAR trt visit baseline y;
RUN;

/* ステップ2：補完データごとに PROC MIXED で解析 */
PROC MIXED DATA=imputed;
  BY _IMPUTATION_;
  CLASS trt visit subject;
  MODEL y = trt visit trt*visit baseline / DDFM=SATTERTH;
  REPEATED visit / SUBJECT=subject TYPE=UN;
  LSMEANS trt*visit / DIFF;
  ODS OUTPUT Diffs=Diffs;
RUN;

/* ステップ3：PROC MIANALYZE でRubin's Rulesを適用 */
PROC MIANALYZE DATA=Diffs;
  BY visit;
  MODELEFFECTS Estimate;
  STDERR StdErr;
RUN;
📝 解釈・補足
FCSステートメントで変数ごとに代入モデルを個別指定できます。MONOTONE REGは単調欠測パターンに特化した高速アルゴリズムです。PROC MIANALYZEのBY visitにより、各時点での治療群間差を統合推定として取り出せます。

モンテカルロシミュレーション：手法比較

完全ケース解析（CCA）・単一代入（平均代入）・多重代入法（MI）のバイアスと標準誤差を比較します。

set.seed(2024)
n_sim  <- 1000
n      <- 100
true_b <- 2.0
results <- matrix(NA, nrow = n_sim, ncol = 3,
                  dimnames = list(NULL, c("CCA","Mean_Imp","MI")))
library(mice)

for (s in 1:n_sim) {
  x   <- rnorm(n)
  y   <- true_b * x + rnorm(n)
  miss_prob <- plogis(-1 + 0.8 * x)
  miss_idx  <- runif(n) < miss_prob
  y_obs     <- y
  y_obs[miss_idx] <- NA
  df <- data.frame(x = x, y = y_obs)

  # CCA
  results[s, "CCA"] <- coef(lm(y ~ x, data = df, na.action = na.omit))["x"]

  # 平均代入
  df_mean <- df
  df_mean$y[miss_idx] <- mean(df$y, na.rm = TRUE)
  results[s, "Mean_Imp"] <- coef(lm(y ~ x, data = df_mean))["x"]

  # MI
  imp <- mice(df, m = 20, printFlag = FALSE, seed = s)
  fit <- with(imp, lm(y ~ x))
  results[s, "MI"] <- summary(pool(fit))[2, "estimate"]
}

bias <- colMeans(results) - true_b
rmse <- sqrt(colMeans((results - true_b)^2))

summary_df <- data.frame(
  Method = c("完全ケース解析","平均代入","多重代入法"),
  Bias   = round(bias, 4),
  RMSE   = round(rmse, 4)
)
print(summary_df)
>       Method    Bias   RMSE
1  完全ケース解析  -0.0031  0.1842
2      平均代入  -0.4127  0.4263
3    多重代入法  -0.0018  0.1798
📝 解釈・補足
MARが成立する条件では、CCAとMIはほぼ不偏（バイアス≈0）ですが、平均代入はバイアス −0.41 と大きく歪んでいます。RMSEはMIが0.1798とCCA（0.1842）よりわずかに小さく、欠測情報を有効活用できていることを示しています。平均代入は推定バイアスが大きく、実務での使用は避けるべきです。

よくある質問（FAQ）

質問回答
代入回数（m）は何回が適切？欠測割合に応じますが、現在の推奨は20〜100回です。欠測率30%なら最低30回を目安にしてください。
解析モデルに含まない変数も代入モデルに入れるべき？はい。補助変数（auxiliary variables）として欠測メカニズムに関連する変数を追加すると推定の精度が向上します。
MNARの場合はMIは使えない？MIはMAR前提ですが、delta-adjustment等の感度分析を組み合わせることでMNAR仮定下での頑健性を評価できます。
代入後の診断はどうやって行う？Rでは plot(imp) でtrace plotを確認し、収束を視覚的に確認します。代入値の分布が観測値と大きくずれていないかも確認してください。
規制提出（CTD）でMIを使う際の注意点は？ICH E9(R1)に従い、主解析の仮定（MAR等）と感度分析の計画をSAPに事前規定することが求められます。代入モデルの詳細も記述します。
SASとRで結果が微妙に異なるのはなぜ？乱数シード・デフォルトの代入アルゴリズムの差異が主な原因です。シードを固定し、アルゴリズム設定を揃えることで一致度が高まります。

実務でのポイント

🔑 まとめ・実務ポイント
  • 受動的代入：派生変数（BMI等）には必ずpassive imputationを設定し、内部整合性を保つ
  • 交互作用項の congenial性：解析モデルの積項は代入モデルにも含める
  • 感度分析の事前計画：MNAR仮定下のdelta-adjustmentをSAPに記載し、規制対応を万全に
  • 補助変数の活用：解析モデル外でも欠測機構に関連する変数は代入モデルに追加する
  • 収束診断：trace plotでmixingを確認し、代入モデルが収束していることを必ず確認する

📚 この記事をより深く理解するための参考書籍

統計・生物統計をさらに深く学びたい方に、おすすめの書籍をご紹介します。

『臨床試験ハンドブック デザインと統計解析』丹後俊郎・上坂浩之（朝倉書店）

製薬企業の生物統計家が押さえるべき臨床試験設計と統計解析を体系的に学べる一冊です。欠測データ処理・感度分析の実務的な観点が充実しており、本記事のMNAR感度分析の背景を理解するうえで最適です。

『臨床試験の事典』丹後俊郎（朝倉書店）

多重代入法・MNAR・ICHガイドラインなど、臨床試験のキーワードを見開き2〜4ページで簡潔に解説しています。辞書的に手元に置いておくと、実務での疑問を素早く解消できます。

関連記事・次のステップ

まとめ

本記事では多重代入法の応用編として、4つのトピックを解説しました。

MICEパッケージの高度な設定では、受動的代入と交互作用項の扱いを取り上げました。派生変数の整合性を保ち、解析モデルとの congenial性を確保することが精度の高い推定につながります。

MNAR感度分析（delta-adjustment）では、MAR仮定が成立しない場合でも、ICH E9(R1)の枠組みに沿って結果の頑健性を評価できることを示しました。

SASのPROC MI（FCS）+ PROC MIANALYZE × PROC MIXEDの連携は、縦断データを扱う製薬業務で直接活用できる実装パターンです。

モンテカルロシミュレーションでは、平均代入のバイアスの大きさと、多重代入法の推定精度の優位性を定量的に確認しました。これらのスキルを組み合わせることで、欠測データを含む臨床試験解析の品質を大幅に向上させることができます。規制提出を視野に入れた感度分析計画も含め、実務の強みにしていただければと思います。

次回は反復測定混合モデル（MMRM）と多重代入法を組み合わせた解析パターンを紹介します。

R SAS 多重代入法 感度分析 欠測データ 臨床試験
ABOUT ME
tomokichi
外資系製薬会社で生物統計家として働ている1児のパパ。生物統計家とは何か、どのようなスキルが必要か、何を行っているのかを共有していきたいと思っております！生物統計に関する最新情報を皆様にお届けすべく、日々奮闘中です。趣味は筋トレ、温泉巡り、家族と散歩。
BLOG：https://www.biostat4869.com