研究者: 3群以上（例：「標準治療 vs 治療A vs 治療B」）ある治療群の生存時間データで交絡調整を行いたいのですが、傾向スコア重み付け（IPTW）のコードや重み付きカプランマイヤー曲線の描き方が分かりません……。
統計ER: 3群以上の多群（Multi-treatment）になっても身構える必要はありません！
WeightIt+cobalt+survival+survminerを使えば、2群比較とほとんど変わらないシンプルな記述で一般化傾向スコアの重み計算からKM曲線描画・Cox回帰（ハザード比算出）まで完結できますよ！さっそく解説しましょう！
臨床研究や観察研究（リアルワールドデータ解析）において、治療法や曝露因子が3群以上（例：「標準治療 vs 新薬A vs 新薬B」「低用量 vs 中用量 vs 高用量」など）存在し、交絡因子を調整した上でカプランマイヤー（KM）生存曲線を描画したりCox比例ハザードモデルでハザード比（HR）を算出したい場面は非常に多い。
2群間の傾向スコア逆確率重み付け（IPTW: Inverse Probability of Treatment Weighting）の手順はWeb上に多く存在するが、3群以上の多群（Multi-treatment）比較になると、傾向スコアの計算や重み付けのコードで途端に詰まってしまう研究者が少なくない。
「多項ロジスティック回帰で一般化傾向スコアを手算出しなければならないのか？」と身構える必要はない。Rの強力なパッケージ群を組み合わせれば、3群以上の生存時間解析であっても、2群比較とほとんど変わらないシンプルな記述で完結できる。
本記事では、WeightIt・cobalt・survival・survminer パッケージを活用し、3群以上の生存時間データにおける重み計算・バランス評価・重み付きKM曲線描画・重み付きCox回帰（HR算出）までを一気通貫で解説する。
3群以上における傾向スコア（一般化傾向スコア：GPS）とIPTWの基本ルール
2群比較における傾向スコアは「治療群に割り当てられる条件付き確率（二項ロジスティック回帰）」であったが、治療法が3群以上（$K \ge 3$）に増えた場合、一般化傾向スコア（Generalized Propensity Score: GPS）の概念が適用される。
GPSでは、多項ロジスティック回帰モデル（Multinomial Logistic Regression）などを用いて、患者の共変量 $X$ のもとで各治療群 $k$ を受ける事後確率 $P(A = k \mid X)$ を算出する。
そして、各症例に対して「実際に受けた治療群 $k$ に割り当てられた確率の逆数」を重み（ATE重み: Inverse Probability Weight）として割り当てる。$$w_i = \frac{1}{P(A = a_i \mid X_i)}$$
数式は複雑に見えるが、Rにおいてはパッケージが裏側で自動的に多項モデルをフィッティングしてくれるため、ユーザーが手動で確率計算を行う必要はない。
使用する神パッケージ4選とパイプラインの役割
3群以上の生存時間IPTW解析をスムーズに行うため、以下の4つのパッケージを連携させる。
WeightIt: 治療変数が3群以上の因子型（factor）であっても、自動的に多項モデルを適用して一括でIPTW重みを計算してくれる万能パッケージ。cobalt: 3群間の全ペア（Pairwise）における標準化平均差（SMD）を一括評価する。survival: 重み付きカプランマイヤー推定量（survfit）および重み付きCox比例ハザードモデル（coxph）を実行する。survminer: 論文掲載レベルの美しく整った多群カプランマイヤー曲線（ggsurvplot）を描画する。
【実践ステップ】Rで3群以上の生存時間データをIPTW解析する全コード
ここからは、実際にR環境で動作するサンプルコードを用いて、Step 1からStep 5まで順番に解説する。
まずは必要なパッケージを読み込んでおこう。
# 必要なパッケージの読み込み
library(tidyverse)
library(WeightIt)
library(cobalt)
library(survival)
library(survminer)
Step 1: サンプルデータ（3群の生存時間データ）の準備
動作確認用として、3つの治療群（Group_A, Group_B, Group_C）と生存時間（time）、イベント有無（status）、共変量を含むサンプルデータを作成する。
# サンプルデータの作成（3群比較）
set.seed(123)
n <- 400
df_surv <- tibble(
id = 1:n,
# 3群の治療変数（Factor型で定義）
treat = factor(sample(c("Group_A", "Group_B", "Group_C"), n, replace = TRUE, prob = c(0.4, 0.3, 0.3))),
age = rnorm(n, mean = 65, sd = 10),
sex = factor(sample(c("Male", "Female"), n, replace = TRUE)),
bmi = rnorm(n, mean = 24, sd = 4),
sbp = rnorm(n, mean = 135, sd = 15),
# 追跡時間とイベント（1: 発生, 0: 打切り）
time = round(runif(n, min = 1, max = 60), 1),
status = sample(c(0, 1), n, replace = TRUE, prob = c(0.6, 0.4))
)
head(df_surv)
Step 2: WeightIt による重み付け（一般化傾向スコア）
WeightIt パッケージの weightit() 関数は、目的変数（治療変数）が3群以上の因子型（factor）になっている場合、自動的に多項モデル（method = "glm"）を選択して各群の一般化傾向スコアとIPTW重みを算出する。
# 3群に対するIPTW重みの計算
W_multi <- weightit(
treat ~ age + sex + bmi + sbp,
data = df_surv,
method = "glm", # 多項ロジスティック回帰が内部で適用される
estimand = "ATE" # ATE重みを指定
)
# 重みオブジェクトの確認
summary(W_multi)
2群比較のコードと全く同じ構文で3群の重み付けが完了する点が、WeightIt を使う最大のメリットである。
Step 3: 3群間の共変量バランス評価（cobalt::bal.tab）
3群以上の解析では、各治療群のペアすべて（Group A vs B、Group A vs C、Group B vs C）における標準化平均差（Pairwise SMD） がバランスの評価対象となる。
cobalt パッケージの bal.tab() を使用すれば、全ペアの最大SMD（Max SMD）や統合結果を一括で確認できる。
# 3群間の共変量バランス（SMD）を一括評価
bal.tab(W_multi, un = TRUE, stats = c("m"), thresholds = c(m = 0.1))
出力結果において、重み付け後（Adj）の最大SMD（Max.Diff.Adj）がすべて 0.1 未満 に収まっていれば、3群間での共変量バランス調整は成功である。
Step 4: 重み付きカプランマイヤー（KM）曲線の描画（survminer）
重み付け調整後の生存確率を推定し、カプランマイヤー曲線を描画する。
survival::survfit() の weights 引数に算出された重み（W_multi$weights）を渡すことで、簡単に重み付きKM推定量を算出できる。
# 重み付きカプランマイヤー曲線の推定
fit_km <- survfit(
Surv(time, status) ~ treat,
data = df_surv,
weights = W_multi$weights
)
# survminer による多群KM曲線の描画
ggsurvplot(
fit_km,
data = df_surv,
pval = TRUE, # ログランク検定P値の表示
risk.table = TRUE, # Risk Tableの表示
legend.title = "Treatment",
legend.labs = c("Group A", "Group B", "Group C"),
palette = c("#E41A1C", "#377EB8", "#4DAF4A"), # カラーパレット
title = "Weighted Kaplan-Meier Survival Curves (3 Groups)",
xlab = "Time (Months)",
ylab = "Overall Survival Probability"
)
調整前の背景バイアスが取り除かれた、美しく信頼性の高い3群比較の生存曲線が出力される。
Step 5: 重み付きCox比例ハザードモデルによるハザード比（HR）算出
各治療群間のハザード比（HR）および95%信頼区間を算出するため、重み付きCox比例ハザードモデルを実行する。
IPTWを適用した疑似集団における回帰分析では、標準誤差が過小評価されるのを防ぐため、堅牢な標準誤差（Robust Standard Error / Sandwich Estimator） を指定する必要がある。coxph() 内で robust = TRUE または cluster = id を指定する。
# 重み付きCox比例ハザードモデルの実行
fit_cox <- coxph(
Surv(time, status) ~ treat,
data = df_surv,
weights = W_multi$weights,
robust = TRUE # 堅牢な標準誤差を算出
)
# 解析結果（ハザード比と95%CI）の表示
summary(fit_cox)
出力結果の解釈:
基準群（デフォルトではアルファベット順の第1水準 Group_A）に対する Group_B および Group_C のハザード比（exp(coef)）とその95%信頼区間・P値が出力される。
3群以上のIPTW解析で気をつけたい注意点
1. 極端な重み（Extreme Weights）への対処
群数が3群、4群と増えるにつれ、一般化傾向スコアの推定値が0に極めて近い症例が発生しやすくなり、極端に大きな重み（Extreme Weights） が生じるリスクが高まる。
もし特定症例の重みが突出して大きい場合は、WeightIt 内のトリミング機能（例: trim(W_multi, at = 0.99)）や、安定化重み（Stabilized Weights）の検討を行うことが推奨される。
2. ハザード比の基準群（Reference）の指定
3群比較では「どの治療群を対照（基準）としてハザード比を算出するか」をあらかじめ決定しておく必要がある。
意図した群（例: 標準治療群）を基準に設定したい場合は、解析前に forcats::fct_relevel() や relevel() を使って因子水準（Factor Levels）の第1水準に指定しておこう。
# Group_A を基準群（Reference）に明示的に設定
df_surv <- df_surv %>%
mutate(treat = fct_relevel(treat, "Group_A"))
【まとめ】コピペで動く一括テンプレートコード
最後に、本記事で紹介した「3群以上の生存時間データに対するIPTW解析」の一式コードをまとめた。自身のデータセット名・変数名に書き換えてコピー＆ペーストで活用してほしい。
# === 3群以上の生存時間データ IPTW解析 一括テンプレート ===
library(tidyverse)
library(WeightIt)
library(cobalt)
library(survival)
library(survminer)
# 1. 基準群の整理（必要に応じて基準群を第1水準に設定）
df <- df_surv %>%
mutate(treat = fct_relevel(treat, "Group_A"))
# 2. IPTW重み（一般化傾向スコア）の計算
W <- weightit(
treat ~ age + sex + bmi + sbp, # 治療変数 ~ 共変量
data = df,
method = "glm",
estimand = "ATE"
)
# 3. 3群間の共変量バランス評価（SMD < 0.1 の確認）
bal.tab(W, un = TRUE, thresholds = c(m = 0.1))
# 4. 重み付きカプランマイヤー曲線の推定と描画
fit_km <- survfit(Surv(time, status) ~ treat, data = df, weights = W$weights)
ggsurvplot(fit_km, data = df, pval = TRUE, risk.table = TRUE)
# 5. 重み付きCox回帰（ハザード比の算出）
fit_cox <- coxph(Surv(time, status) ~ treat, data = df, weights = W$weights, robust = TRUE)
summary(fit_cox)
難解に見える多群の生存時間因果推論も、WeightIt ＋ survival ＋ survminer のパイプラインを使えば非常にスマートに処理できる。
適切な交絡調整を行い、説得力の高い臨床研究結果を発表していこう。
コメント