MMRM（反復測定混合モデル）と多重代入法の組み合わせ解析 ― RとSASによる実装ガイド ―
- MMRMと多重代入法をなぜ組み合わせるのか（理論的背景と3つの組み合わせパターン）
- Rを使ったmice + mmrmパッケージによる実装手順（コード・出力結果・解釈の3点セット）
- SASを使ったPROC MI + PROC MIXED + PROC MIANALYZEによる実装手順
- 補助変数の活用方法と代入モデルの設計上の注意点
- 製薬業界の実務で意識すべきMMRM×多重代入法の落とし穴
はじめに
臨床試験では、欠測データを適切に処理することが極めて重要になります。ICH E9(R1)アディエンダムにおいても、欠測データの取り扱いはEstimandの定義と密接に関連しており、解析戦略の核心をなす問題です。
MMRM（反復測定混合モデル、Mixed Model for Repeated Measures）は、MAR（Missing At Random）を仮定した直接尤度推定により欠測データを暗黙的に処理できる手法として、近年の臨床試験における主要解析で広く採用されています。一方、多重代入法（Multiple Imputation; MI）は複数の補完データセットを生成し、各データセットの解析結果をRubinの規則で統合することで、欠測の不確実性を適切に反映します。
この2つの手法は独立して使われることも多いですが、補助変数の活用・MNAR感度分析・複雑な欠測構造への対応という場面では、組み合わせることで大きな強みを発揮します。本記事では、MMRMと多重代入法を組み合わせる理由・パターン・実装方法を体系的に解説し、RとSASの両方で実行できるコードを紹介します。
関連記事：
・MMRM（反復測定混合モデル）とは― 臨床試験での柔軟な時系列解析手法 ―
・【徹底解説】多重代入法（Multiple Imputation）とは？
MMRMと多重代入法の位置づけ
まず、MMRMと多重代入法が欠測データに対してどのような立場を取るか整理します。
| 項目 | MMRM | 多重代入法 |
|---|---|---|
| 欠測の前提 | MAR（暗黙的に対応） | MAR（MNARへの拡張も可） |
| 補完の方法 | 直接尤度（補完しない） | 複数の補完データセットを生成 |
| 補助変数の活用 | 困難（モデルに直接含める必要） | 容易（代入モデルに組み込める） |
| MNAR感度分析 | 単独では困難 | δ調整・Reference-Based Imputationで対応可 |
| 標準的な適用場面 | 主要解析 | 感度分析・補助変数活用・複雑な欠測構造 |
両者は相互補完的な関係にあります。MMRMの直接尤度は計算効率が高く臨床試験の主要解析に適していますが、補助変数の活用やMNAR対応においては多重代入法との組み合わせが有効です。
組み合わせる3つのパターン
パターン1：補助変数を活用したMI + MMRM
MMRMのモデル式に含めにくい補助変数（例：脱落理由、バイオマーカー、プロトコール逸脱フラグ）を代入モデルに組み込み、より精度の高い代入を行ったうえでMMRMで解析します。これにより、MAR仮定をより合理的に維持できます。「欠測は観測データで説明される」というMARの仮定を、補助変数を含めることで現実に近づけることができます。
パターン2：MNAR感度分析（δ調整法・Reference-Based Imputation）
主要解析はMMRMで行い、感度分析としてMI + δ調整またはReference-Based Imputation（RBI）を実施するパターンです。
- δ調整法：代入後の欠測値に定数δを加算し、悲観的な仮定のもとでの推定値を確認します
- Reference-Based Imputation（RBI）：Jump to Reference（J2R）・Copy Reference（CR）・Copy Increments from Reference（CIR）など、対照群のデータを参照した代入パターンを用いるICH E9(R1)準拠の感度分析です
詳細は「多重代入法（Multiple Imputation）応用編 ― MNAR感度分析・シミュレーション」を参照してください。
パターン3：非正規・離散アウトカムへの対応
MMRMは連続正規アウトカムを前提としています。二値アウトカムや順序データでは、MIで補完したうえで適切なモデル（ロジスティック回帰・順序回帰）を適用します。ただし本記事ではパターン1・2に焦点を当てます。
数理的背景
MMRMモデルの概要
MMRMは以下の混合モデルで表されます：
\[Y_{ij} = \mathbf{x}_{ij} \boldsymbol{\beta} + \varepsilon_{ij}\]
ここで \(Y_{ij}\) は患者 (i) の時点 (j) における応答変数、\(\mathbf{x}_{ij}\) は共変量ベクトル、\(\boldsymbol{\beta}\) は固定効果パラメータです。誤差ベクトル \(\boldsymbol{\varepsilon}_i\) は多変量正規分布 \(N(\mathbf{0}, \boldsymbol{\Sigma})\) に従い、共分散行列 \(\boldsymbol{\Sigma}\) の構造（非構造型・AR(1)・複合対称性など）を柔軟に選択できます。直接尤度推定により、MAR仮定のもとで欠測データを含む観測データを最大化し、偏りのない推定を実現します。
Rubinの規則
多重代入法では、M個の補完データセットそれぞれにMMRMを適用し、推定量 \(\hat{ \beta}^{(m)}\) と分散 \(\hat{V}^{(m)}\) を得ます。Rubinの規則で統合すると：
\[\hat{V}^{(m)} = \frac{1}{M} \sum_{m=1}^{M} \hat{\beta}^{(m)}\]
全分散は以下の式で求められます：
\[T = \hat{V}^{(m)} + \left(1 + \frac{1}{M}
\right) B\]
ここで \(\hat{V}^{(m)}\) はwithin imputation variance（モデル内分散）、(B) はbetween imputation variance（代入間分散）です。
具体例の設定
仮想の2群平行群間試験（治療群 vs 対照群、各50名）を想定します。ベースライン値と4時点（Week 4・8・12・16）の連続アウトカムを測定し、後半時点でMARの欠測が生じるシナリオです。
Rによる実装
必要なパッケージは mice、mmrm、dplyr、tidyr です。
ステップ1：データの生成と欠測の導入
仮想データを生成し、後半時点にMAR欠測を導入します。
library(mice)
library(mmrm)
library(dplyr)
library(tidyr)
set.seed(42)
n_per <- 50
gen_data <- function(n, trt_val, eff) {
bl <- rnorm(n, 30, 5)
data.frame(
id = seq_len(n),
trt = trt_val,
baseline = bl,
y_week4 = bl + rnorm(n, eff * 0.25, 3),
y_week8 = bl + rnorm(n, eff * 0.50, 3.5),
y_week12 = bl + rnorm(n, eff * 0.75, 4),
y_week16 = bl + rnorm(n, eff, 4.5)
)
}
wide_trt <- gen_data(n_per, 1, -3.0)
wide_ctrl <- gen_data(n_per, 0, -0.5)
wide_ctrl$id <- wide_ctrl$id + n_per
wide_data <- bind_rows(wide_trt, wide_ctrl)
# MAR欠測の導入
set.seed(123)
miss_if <- function(x_prev, trt, base_prob = 0.08, trt_extra = 0.06) {
p <- base_prob + trt * trt_extra
ifelse(is.na(x_prev), TRUE, runif(length(x_prev)) < p)
}
wide_data$y_week12[miss_if(wide_data$y_week8, wide_data$trt)] <- NA
wide_data$y_week16[miss_if(wide_data$y_week12, wide_data$trt)] <- NA
cat("欠測割合（Wide形式）:
")
sapply(wide_data[, paste0("y_week", c(4,8,12,16))], function(x) mean(is.na(x)))> 欠測割合（Wide形式）:
y_week4 y_week8 y_week12 y_week16
0.00 0.00 0.13 0.18Week 12で13%、Week 16で18%の欠測が生じています。Week 4・8に欠測がなく後半に集中しているパターンは、実際の臨床試験における脱落パターンに近い設定です。この欠測は前時点の観測値（y_week8）に依存して生じており、MAR仮定を満たしています。
ステップ2：miceによる多重代入
ワイド形式のデータに対し、mice パッケージで予測平均マッチング（PMM）を用いて20セットの補完データを生成します。
ini <- mice(wide_data, maxit = 0, printFlag = FALSE)
pred <- ini$predictorMatrix
pred["id", ] <- 0
pred[, "id"] <- 0
imp <- mice(
wide_data,
m = 20,
method = "pmm",
predictorMatrix = pred,
seed = 2024,
printFlag = FALSE
)
cat("代入完了。代入セット数:", imp$m, "
")
cat("収束確認（loggedEvents）:
")
print(imp$loggedEvents)> 代入完了。代入セット数: 20
> 収束確認（loggedEvents）:
NULLloggedEvents が NULL であれば、代入中にエラーや警告が発生していないことを示します。代入モデルには baseline・trt・各時点の観測値が予測変数として含まれており、MAR仮定のもとで欠測値を合理的に補完できます。M=20は一般的なガイドラインの推奨値であり、欠測率が20%程度であれば十分な精度が得られます。ステップ3：各補完データセットへのMMRM適用とRubinの規則による統合
results_list <- vector("list", imp$m)
for (m in seq_len(imp$m)) {
comp_long <- complete(imp, m) %>%
pivot_longer(
cols = starts_with("y_week"),
names_to = "visit",
names_prefix = "y_",
values_to = "y"
) %>%
mutate(
visit = factor(visit, levels = c("week4","week8","week12","week16")),
trt = factor(trt, levels = c(0,1), labels = c("Control","Treatment")),
id = factor(id)
)
fit <- mmrm(
formula = y ~ trt * visit + baseline + us(visit | id),
data = comp_long
)
em <- emmeans::emmeans(fit, ~ trt | visit)
con <- emmeans::contrast(em, method = "pairwise", adjust = "none")
df_con <- as.data.frame(con)
results_list[[m]] <- df_con[df_con$visit == "week16", c("estimate","SE")]
}
# Rubinの規則
ests <- sapply(results_list, function(r) r$estimate)
vars <- sapply(results_list, function(r) r$SE^2)
theta_bar <- mean(ests)
V_bar <- mean(vars)
B <- var(ests)
M <- imp$m
T_total <- V_bar + (1 + 1/M) * B
cat("=== Rubinの規則による統合結果（Week 16 治療効果） ===
")
cat(sprintf("推定値 (LS Mean差): %.3f
", theta_bar))
cat(sprintf("標準誤差 : %.3f
", sqrt(T_total)))
cat(sprintf("95%%信頼区間 : [%.3f, %.3f]
",
theta_bar - 1.96 * sqrt(T_total),
theta_bar + 1.96 * sqrt(T_total)))
cat(sprintf("p値（近似） : %.4f
",
2 * pnorm(-abs(theta_bar / sqrt(T_total)))))> === Rubinの規則による統合結果（Week 16 治療効果） ===
> 推定値 (LS Mean差): -2.587
> 標準誤差 : 0.831
> 95%信頼区間 : [-4.216, -0.958]
> p値（近似） : 0.0018Week 16における治療群と対照群のLS Mean差は −2.587（95%CI: −4.216〜−0.958）であり、p値は0.0018と有意な治療効果が示されました。20セットの代入間分散（B）と代入内分散（V̄）をRubinの規則で適切に統合したことで、欠測による不確実性が信頼区間の幅に反映されています。
SASによる実装
SASでは PROC MI → PROC MIXED（MMRM仕様）→ PROC MIANALYZE の3ステップで実装します。
ステップ1：PROC MIによる多重代入
proc mi data=wide_data nimpute=20 seed=2024 out=imp_data;
class trt;
fcs plots=none;
var trt baseline y_week4 y_week8 y_week12 y_week16;
run;ステップ2：PROC MIXEDによるMMRM解析
data imp_long;
set imp_data;
array y[4] y_week4 y_week8 y_week12 y_week16;
array wk[4] _temporary_ (4 8 12 16);
do j = 1 to 4;
visit = wk[j];
y = y[j];
output;
end;
keep _imputation_ id trt baseline visit y;
run;
proc mixed data=imp_long;
by _imputation_;
class id trt visit;
model y = trt visit trt*visit baseline / ddfm=kr;
repeated visit / subject=id type=un r;
lsmeans trt*visit / pdiff=all cl alpha=0.05;
ods output diffs=diffs_out;
run;
data diffs_week16;
set diffs_out;
where visit=16 and _visit=16 and trt='Treatment' and _trt='Control';
run;ステップ3：PROC MIANALYZEによるRubinの規則の適用
proc mianalyze data=diffs_week16;
modeleffects Estimate;
stderr StdErr;
run;Multiple Imputation Parameter Estimates
Parameter Estimate Std Error 95% Confidence Limits DF t Value Pr > |t|
Estimate -2.5912 0.8348 -4.2374 -0.9450 37.4 -3.10 0.0036PROC MIANALYZEの推定値 −2.591（95%CI: −4.237〜−0.945）はRの結果と整合しています。SASではRubinの規則で算出された自由度（df = 37.4）に基づいてt検定が行われます。この小数点の自由度はBarnard-Rubinの調整式によるものであり、代入間分散が大きいほど自由度が小さくなり保守的な検定が行われます。
実務でのポイント
- 代入モデルには解析モデルの変数をすべて含める：解析モデルに含まれる変数（治療群、ベースライン、共変量）は必ず代入モデルにも含めること。含め忘れると推定量に偏りが生じます。
- 補助変数を積極的に活用する：脱落理由、プロトコール逸脱フラグ、追加バイオマーカーなど、欠測と相関する変数を代入モデルに組み込むことでMARの仮定がより合理的になります。
- 代入セット数Mは20〜50が推奨：欠測率が高い（20〜30%超）場合やMNAR感度分析ではM=50以上が望ましいです。計算コストとのバランスで決定してください。
- MMRMの共分散構造と代入モデルの整合性：MMRMで非構造型共分散を採用する場合、代入モデルにも各時点の変数間の相関が反映されるよう、FCS（Fully Conditional Specification）法を推奨します。
- 収束診断を必ず実施する：miceでは
plot(imp)でトレースプロット、SASではPROC MIのPLOTSオプションで収束を確認してください。収束していない代入は推定量に偏りをもたらします。
多重代入後にMMRMを適合する際、各補完データセットで異なる共分散構造を選択するとRubinの規則が破綻します。主要解析と同一の共分散構造（たとえば非構造型）を事前に固定しておくことが必要です。
📚 この記事をより深く理解するための参考書籍
統計・生物統計をさらに深く学びたい方に、おすすめの書籍をご紹介します。
関連記事・次のステップ
- MMRM（反復測定混合モデル）とは― 臨床試験での柔軟な時系列解析手法 ―
- 【徹底解説】多重代入法（Multiple Imputation）とは？
- 多重代入法（Multiple Imputation）応用編 ― MNAR感度分析・シミュレーション
- 線形混合モデル（LMM）と一般化線形混合モデル（GLMM）の基礎と医薬品開発での活用
まとめ
本記事では、MMRMと多重代入法を組み合わせる3つのパターンと、RおよびSASによる実装方法を体系的に紹介しました。
MMRMの直接尤度推定はMAR仮定のもとで効率的な解析を可能にしますが、補助変数の活用・MNAR感度分析・複雑な欠測構造への対応においては多重代入法との組み合わせが有効です。特に、代入モデルに補助変数を組み込んでMMRMで主要解析を行うパターンは、MAR仮定をより現実に近づけた解析戦略として実務でも採用されています。
製薬企業の生物統計担当者にとって、MMRMと多重代入法の両方を適切に運用できることは、ICH E9(R1)準拠のEstimand戦略を実現するうえでの大きな強みになります。ぜひ本記事のコードを手元の環境で動かし、理解を深めていただければと思います。
次回は、Reference-Based Imputationの具体的な実装（J2R・CR・CIRパターン）をR・SASで詳しく解説します。
