RでIPTW法、Overlap weighting法を用いたKaplan-meier曲線を作成する

こんにちは、eURO2023です！

趣味（？）でRstudioをしています。

Kaplan-Meyer曲線について、前回お話しました。

応用編として、レトロスペクティブなデータ解析で有用なIPTW法を用いる機会があったので、備忘録としてまとめておきます。

更に、最近Overlap weighting法という重み付け法もあるようなので、合わせてまとめてみました！

医師が統計を学ぶ上で、絶対読んだほうがよい！と思われる書籍を上げておきます。

特に、新谷先生が執筆された「今日から使える医療統計」はおすすめです！

私は第一版を持っていますが、第二版が登場したようですね。

とても読みやすく、医療統計を始める上でも、論文の中身を理解する上でも非常に参考になります。

ある程度理解できるようになった今でも、よく読み返す本です。

「医療統計の基礎のキソ」は、更にわかりやすく、

「p値とはなにか？」

という超基本的な内容から、深く掘り下げてわかりやすく伝えてくれています！

第1巻から3巻までありますが、すぐに読めますよ。

下準備

今回も、前回作成した架空のデータシート”X”を用いていきます。

詳細は、下記ブログを参考にしてください。

RでKaplan-Meier曲線を描く方法

Rstudioを用いた、カプランマイヤー曲線の作り方を解説。

urology-tweet.com

2025.02.23

また、既存のCSVファイルやExcelファイルを読み込む場合は、readrパッケージやopenxlsxパッケージを用います。

一般的には、Rの解析にはExcelよりCSVの方が柔軟性があって良いらしいです。が、Excelファイルでもそんなに困ることは（私の場合）ありません。

# CSVファイルの読み込み
install.packages("readr")
library(readr)
X <- read_csv("絶対パス")

# Excelファイルの読み込み
install.packages("openxlsx")
library(openxlsx)
X <- read.xlsx("ファイルの絶対パス")

応用編①：IPTW法によるKaplan-Meier曲線の作成

IPTW法は、”Inverse Probability of Treatment Weighting”（逆確率重み付け法）法の略です。

観察研究において、治療群と対照群のバイアスを調整し、ランダム化試験に近い比較を行うための方法です。

観察研究では、治療を受けるかどうかはランダムに選ばれることはありません。

例えば、高齢者の方が特定の治療を受けやすい、重症患者がより積極的な治療を受ける、など様々なバイアスが実臨床では出てきます。

そのまま治療群と非治療群を比較すると、このようなバイアス、または交絡因子（年齢、性別、BMI など）が結果に影響してしまいます。観察研究が、ランダム化比較試験に比べて研究の質が落ちるとされているのは、このバイアスによる面が大きいと言われています。

IPTW 法を使うことで、治療の割り付けがランダム化されたかのようにデータを調整できると言われており、近年よく用いられるようになってきました。

さて、IPTW法を用いてKaplan-Meier曲線を作成するためのRスクリプトを確認しましょう。

# Propensity Score (PS) モデルの構築
ps_model <- glm(drug_A ~ age + sex + BMI + Hb + PS, data = X, family = binomial)

# PSの推定値を取得
X$ps <- predict(ps_model, type = "response")

# IPTW法の重み計算
X$iptw_weight <- ifelse(X$drug_A == 1, 1 / X$ps, 1 / (1 - X$ps))

# サバイバルオブジェクトの作成
surv_obj <- Surv(time = X$survival_month, event = X$survival)

# IPTW法を適用した生存分析モデル (survival model) の構築
km_fit_iptw <- survfit(Surv(survival_month, survival) ~ drug_A, weights = iptw_weight, data = X)

# IPTW法によるKaplan-Meier曲線の描画
ggsurvplot(
  km_fit_iptw, 
  conf.int = TRUE,               # 信頼区間の表示
  pval = TRUE,                   # Log-rank検定のp値
  risk.table = TRUE,             # リスクテーブル
  title = "Kaplan-Meier Curve (IPTW)",
  legend.labs = c("No Drug A", "Drug A"), 
  legend.title = "Treatment",
  xlab = "Time (months)",
  ylab = "Survival Probability",
  palette = c("#E7B800", "#2E9FDF")
)

Propensity Score (PS) モデルの構築

# Propensity Score (PS) モデルの構築
ps_model <- glm(drug_A ~ age + sex + BMI + Hb + PS, data = X, family = binomial)

glm() は 一般化線形モデル（Generalized Linear Model, GLM） を作成する関数です。

glm()は、Rの基本的なパッケージに含まれているので、インストールやlibrary()の必要はありません！

上記ステータスを入力し、新たなロジスティック回帰モデル“ps_model”を作成します。

ここで一つ注意点があります。説明変数には今回、age, sex, BMI, Hb, PS (performance status)の5項目を入力しました。

この説明変数には明確な数の制限はありません。なので、必要と思われる説明変数は全て入れることができます。

ただ、IPTW法では下記に示すように、外れ値があるとps値が極端に大きくなってしまいます。

説明変数が多くなると、この極端な外れ値が生じるリスクも高まってしまいます。

入れる説明変数については十分に検討が必要です。

# PSの推定値を取得
X$ps <- predict(ps_model, type = "response"

predict(ps_model, type = “response”) は、作成したロジスティック回帰モデル ps_model を使って、各患者の drug_A を受ける確率（傾向スコア）を予測します。

type = “response” を指定することで、確率（0〜1の値） として出力されます。
X$ps <-：データフレーム X に新しい列 ps を追加し、各患者の傾向スコアを格納します。

これによって、データフレームXにpropensity score（どの患者が drug_A を受ける可能性が高いか）を追加することができました！

IPTW法の重み計算

# IPTW法の重み計算
X$iptw_weight <- ifelse(X$drug_A == 1, 1 / X$ps, 1 / (1 - X$ps))

IPTW では、各個人に対して 「治療を受ける確率（Propensity Score; PS）」の逆数を重みとして与えます。

この調整によって、治療を受ける可能性が高い人と低い人の影響をバランスよく調整できます。

X$iptw_weight → データ X に新しい列 iptw_weight を追加。

ifelse(条件, 値1, 値2) → 条件に応じて異なる値を計算します。

X$drug_A == 1（もしdrug_Aが”1″、つまり治療群の場合）→ 1 / X$ps

X$drug_A == 0（もしdrug_Aが”0″、つまり非治療群の場合）→ 1 / (1 – X$ps)

このようにIPTWの重みを計算して、データフレームXに新たに作成したiptw_weightに保存されます。

このように、psまたは1-psの逆数を重み付けとします。

psは0から1の値を取ります。つまりpsが0か1に近くなるにつれて、psは極端に大きくなってしまうので、IPTW法においては特に注意が必要です！

ちなみに、このifelse関数は本当によく使います！Google検索すると沢山解説がありますので、是非確認して使いこなしていただければと思います。

サバイバルオブジェクトの作成

# サバイバルオブジェクトの作成
surv_obj <- Surv(time = X$survival_month, event = X$survival)

このコマンドで、生存オブジェクト（Survival Object） を作成します。

この部分は、前回ブログにも出てきましたね。

詳しくは前回のブログにも記載しているので、確認してみてください。

RでKaplan-Meier曲線を描く方法

Rstudioを用いた、カプランマイヤー曲線の作り方を解説。

urology-tweet.com

2025.02.23

IPTW法を適用した生存分析モデル (survival model) の構築

# IPTW法を適用した生存分析モデル (survival model) の構築
km_fit_iptw <- survfit(Surv(survival_month, survival) ~ drug_A, weights = iptw_weight, data = X)

このコマンドで、IPTW法を適用したKaplan-Meier 生存曲線 (Kaplan-Meier Survival Curve) を作成するための 生存分析モデル (survival model) を構築します。

基本的には、IPTW法を用いない、生データでのkm_fitを作成する場合と変わりません（前回ブログを参照してください）。

ここにweightsを加えることで、IPTW法を適用できます！

IPTW法によるKaplan-Meier曲線の描画

# IPTW法によるKaplan-Meier曲線の描画
ggsurvplot(
  km_fit_iptw, 
  conf.int = TRUE,               # 信頼区間の表示
  pval = TRUE,                   # Log-rank検定のp値
  risk.table = TRUE,             # リスクテーブル
  title = "Kaplan-Meier Curve (IPTW)",
  legend.labs = c("No Drug A", "Drug A"), 
  legend.title = "Treatment",
  xlab = "Time (months)",
  ylab = "Survival Probability",
  palette = c("#E7B800", "#2E9FDF")
)

ggsurvplotを用いて、Kaplan-Meier曲線を描きます。

ここも、IPTW法を用いない場合とほぼ一緒です。km_fitをkm_fit_iptwに置き換えるだけです。

data = Xはあってもなくてもどちらでも作図できます。

IPTW法であることをわかりやすくするために、title = “Kaplan-Meier Curve (IPTW)”にてタイトルをつけました。

実際の人数（計100名）と異なり、今回はだいぶ人数が増えた状態で解析されています。

これは、IPTW法を適用した場合、実際のサンプルサイズ（n数）とは異なる “有効サンプルサイズ（effective sample size, ESS）” で解析されるからのようです。

risk.table = FALSEとすれば、リスクテーブルが表示されなくなるので、ややっこしくなくて良いかもしれませんね！

ChatGPTに尋ねたところ、IPTW法には以下のメリット・デメリットがあるようです。

✅ メリット

• 交絡因子を調整し、ランダム化試験に近い比較ができる
• すべての共変量のバランスを同時に取れる（マッチング法より柔軟）

❌ デメリット

• 極端なPS（0 に近い or 1 に近い）があると、重みが大きくなりすぎる
• 適切なPS モデルを選ばないと、かえってバイアスが残る

応用編②：Overlap weighting法によるKaplan-Meier曲線の作成

Overlap weighting法は、IPTW法の欠点を補完する重み付け推定法として、最近注目されています。

IPTW法と同様、傾向スコア（Propensity Score, ps）を用いた重み付け推定の一つです。

治療効果を推定する際に、処置群（治療を受けた群）と対照群（治療を受けなかった群）の比較をより公平に行うための方法です。

特に、「オーバーラップ（重なり）」が大きい部分に重点を置くことで、より信頼性の高い推定を可能にします。従来のIPTW（Inverse Probability of Treatment Weighting）法と比べて、極端な重みがつきにくく、推定の安定性が向上する点が特徴です。

以下のブログに詳しく解説がなされており、参考にさせていただきました！

library(survey)
library(survival)
library(survminer)

# Propensity Score (PS) モデルの構築
ps_model <- glm(drug_A ~ age + sex + BMI + Hb + PS, data = X, family = binomial)

# PSの推定値を取得
X$ps <- predict(ps_model, type = "response")

# Overlap Weighting法の重み計算
X$ow_weight <- ifelse(X$drug_A == 1, X$ps, 1 - X$ps)

# サバイバルオブジェクトの作成
surv_obj <- Surv(time = X$survival_month, event = X$survival)

# Overlap weighting法を適用した生存分析モデル (survival model) の構築
km_fit_ow <- survfit(Surv(survival_month, survival) ~ drug_A, weights = ow_weight, data = X)

# Overlap weighting法によるKaplan-Meier 曲線の描画
ggsurvplot(
  km_fit_ow, 
  conf.int = TRUE,
  pval = TRUE,
  risk.table = TRUE,
  title = "Kaplan-Meier Curve (Overlap Weighting)",
  legend.labs = c("No Drug A", "Drug A"), 
  legend.title = "Treatment",
  xlab = "Time (months)",
  ylab = "Survival Probability",
  palette = c("#FC4E07", "#00AFBB")
)

Propensity Score (PS) モデルの構築

# Propensity Score (PS) モデルの構築
ps_model <- glm(drug_A ~ age + sex + BMI + Hb + PS, data = X, family = binomial)

# PSの推定値を取得
X$ps <- predict(ps_model, type = "response")

ここの部分はIPTW法と全く同じですね。

Overlap Weighting法の重み計算

# Overlap Weighting法の重み計算
X$ow_weight <- ifelse(X$drug_A == 1, X$ps, 1 - X$ps)

Overlap Weighting法では、各個体に以下のような重み（weight）を付けます：

IPTWでは、PSが0または1に近いと極端に大きな重みがつくことがあります（不安定）。

Overlap Weightingでは、PSが0または1に近い個体の重みが小さくなるため、安定した推定が可能です。

PSが0.2～0.8程度の中間層に重点を置くことで、バランスが取りやすくなると言われています。

サバイバルオブジェクトの作成

# サバイバルオブジェクトの作成
surv_obj <- Surv(time = X$survival_month, event = X$survival)

この部分も、IPTW法と同じです。

Overlap weighting法を適用した生存分析モデル (survival model) の構築

# Overlap weighting法を適用した生存分析モデル (survival model) の構築
km_fit_ow <- survfit(Surv(survival_month, survival) ~ drug_A, weights = ow_weight, data = X)

ここも、weightをoverlap weighting法で得た重み付け値にする以外は、IPTW法と同じです。

Overlap weighting法によるKaplan-Meier曲線の描画

# Overlap weighting法によるKaplan-Meier曲線の描画
ggsurvplot(
  km_fit_ow, 
  conf.int = TRUE,
  pval = TRUE,
  risk.table = TRUE,
  title = "Kaplan-Meier Curve (Overlap Weighting)",
  legend.labs = c("No Drug A", "Drug A"), 
  legend.title = "Treatment",
  xlab = "Time (months)",
  ylab = "Survival Probability",
  palette = c("#FC4E07", "#00AFBB")
)

ggsurvplotを用いて、Kaplan-Meier曲線を描きます。

ここも、IPTW法とほぼ一緒です。km_fit_iptwをkm_fit_owに置き換えるだけです。

Overlap Weighting法を適用した場合も、実際のサンプルサイズ（n数）とは異なる “有効サンプルサイズ（effective sample size, ESS）” で解析されるため、実際のサンプルサイズとは異なっていることに注意してください。

リスクテーブルを消去したい場合は、risk.table = FALSEとしましょう。

このように、上記の流れで解析をすれば、通常のKaplan-Meier法だけでなくIPTW法も、Overlap Weighting法も、解析対象に応じて柔軟に変更することが可能で便利です！

COX比例ハザードモデル

蛇足ですが、この架空データベースを用いてCox比例ハザードモデルもIPTW法・

# IPTW法でのCox比例ハザードモデル
cox_iptw <- coxph(Surv(survival_month, survival) ~ drug_A + age + sex + BMI + Hb + PS, weights = X$iptw_weight, data = X)
summary(cox_iptw) #結果の表示

# Overlap weighting法でのCox比例ハザードモデル
cox_ow <- coxph(Surv(survival_month, survival) ~ drug_A + age + sex + BMI + Hb + PS, weights = X$ow_weight, data = X)
summary(cox_ow) #結果の表示

例えばIPTW法におけるCOX比例ハザードモデルの結果は、下記のように表示されます。

> summary(cox_iptw) #結果の表示
Call:
coxph(formula = Surv(survival_month, survival) ~ drug_A + age + 
    sex + BMI + Hb + PS, data = X, weights = X$iptw_weight)

  n= 100, number of events= 54 

            coef exp(coef)  se(coef) robust se      z Pr(>|z|)  
drug_A  0.055010  1.056552  0.201710  0.290078  0.190   0.8496  
age     0.018510  1.018682  0.007917  0.011035  1.677   0.0935 .
sex    -0.427481  0.652150  0.210346  0.307120 -1.392   0.1640  
BMI     0.008725  1.008763  0.019872  0.032386  0.269   0.7876  
Hb     -0.054681  0.946787  0.044844  0.062609 -0.873   0.3825  
PS     -0.033024  0.967515  0.073791  0.100393 -0.329   0.7422  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

       exp(coef) exp(-coef) lower .95 upper .95
drug_A    1.0566     0.9465    0.5984     1.866
age       1.0187     0.9817    0.9969     1.041
sex       0.6521     1.5334    0.3572     1.191
BMI       1.0088     0.9913    0.9467     1.075
Hb        0.9468     1.0562    0.8375     1.070
PS        0.9675     1.0336    0.7947     1.178

Concordance= 0.563  (se = 0.044 )
Likelihood ratio test= 11.62  on 6 df,   p=0.07
Wald test            = 8.69  on 6 df,   p=0.2
Score (logrank) test = 11.46  on 6 df,   p=0.08,   Robust = 6.99  p=0.3

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).

さて、Cox比例ハザードモデルにおいて、そもそもIPTW法やOverlap Weighting法を適用することに、一体何の意味があるのでしょうか？

実際にGoogle検索してみると、ちらほら「IPTW Cox回帰」について解説されているサイトをみかけました。

しかし、その意義についてはあまり詳述されていません…。

ということで、ChatGPTさんに聞いてみました（笑）

📌 Cox比例ハザードモデルにIPTW法やOverlap Weighting法を適用する意義

Cox比例ハザードモデルによる多変量解析 では、治療（例: drug_A）と生存期間の関係を評価する際、交絡因子（年齢・性別・BMI など）の影響を調整する必要があります。
この調整方法として、IPTW（Inverse Probability of Treatment Weighting）法や Overlap Weighting（OW）法を適用することが可能です。

---

1️⃣ IPTW法を適用する意義

IPTW法では、各患者の Propensity Score（PS; 治療を受ける確率）を用いて重み付けを行い、疑似的に無作為化比較試験（RCT）に近づけることを目指します。

✅ 交絡因子を調整し、治療効果のバイアスを減らす
✅ 「治療を受ける確率が高い患者」と「治療を受ける確率が低い患者」のバランスを取る
✅ RCTのような「治療の割り付けがランダムな状況」に近づける

2️⃣ Overlap Weighting（OW）法を適用する意義

Overlap Weighting法は、Propensity Scoreが0.5に近い患者（治療と非治療のどちらにも属する可能性が高い患者）に重点を置いて解析する方法 です。

✅ 治療群と非治療群の「共通サポート領域（overlap region）」を重視する
✅ 極端な PS を持つ患者（0 or 1に近い患者）を減らし、外れ値の影響を抑える
✅ バランスの取れた母集団における治療効果を推定できる

とのこと。

意義はありそうですね。特にIPTW法を用いることで、後ろ向き試験であってもRCTのような均一な患者背景に近づけた上で、Cox比例ハザードモデルを適用できる、というメリットがありそうです。

さいごに

# スクリプト全景

library(survey)
library(survival)
library(survminer)

# CSVファイルの読み込み
install.packages("readr")
library(readr)
X <- read_csv("絶対パス")

# Excelファイルの読み込み
install.packages("openxlsx")
library(openxlsx)
X <- read.xlsx("ファイルの絶対パス")

# Propensity Score (PS) モデルの構築
ps_model <- glm(drug_A ~ age + sex + BMI + Hb + PS, data = X, family = binomial)

# PSの推定値を取得
X$ps <- predict(ps_model, type = "response")

# IPTW法の重み計算
X$iptw_weight <- ifelse(X$drug_A == 1, 1 / X$ps, 1 / (1 - X$ps))

# Overlap Weighting法の重み計算
X$ow_weight <- ifelse(X$drug_A == 1, X$ps, 1 - X$ps)

# サバイバルオブジェクトの作成
surv_obj <- Surv(time = X$survival_month, event = X$survival)

# IPTW法を適用した生存分析モデル (survival model) の構築
km_fit_iptw <- survfit(Surv(survival_month, survival) ~ drug_A, weights = iptw_weight, data = X)

# Overlap weighting法を適用した生存分析モデル (survival model) の構築
km_fit_ow <- survfit(Surv(survival_month, survival) ~ drug_A, weights = ow_weight, data = X)

# IPTW法によるKaplan-Meier曲線の描画
ggsurvplot(
  km_fit_iptw, 
  conf.int = TRUE,               # 信頼区間の表示
  pval = TRUE,                   # Log-rank検定のp値
  risk.table = TRUE,             # リスクテーブル
  title = "Kaplan-Meier Curve (IPTW)",
  legend.labs = c("No Drug A", "Drug A"), 
  legend.title = "Treatment",
  xlab = "Time (months)",
  ylab = "Survival Probability",
  palette = c("#E7B800", "#2E9FDF")
)

# Overlap weighting法によるKaplan-Meier曲線の描画
ggsurvplot(
  km_fit_ow, 
  conf.int = TRUE,
  pval = TRUE,
  risk.table = TRUE,
  title = "Kaplan-Meier Curve (Overlap Weighting)",
  legend.labs = c("No Drug A", "Drug A"), 
  legend.title = "Treatment",
  xlab = "Time (months)",
  ylab = "Survival Probability",
  palette = c("#FC4E07", "#00AFBB")
)

# IPTW法でのCox比例ハザードモデル
cox_iptw <- coxph(Surv(survival_month, survival) ~ drug_A + age + sex + BMI + Hb + PS, weights = X$iptw_weight, data = X)
summary(cox_iptw) #結果の表示

# Overlap weighting法でのCox比例ハザードモデル
cox_ow <- coxph(Surv(survival_month, survival) ~ drug_A + age + sex + BMI + Hb + PS, weights = X$ow_weight, data = X)
summary(cox_ow) #結果の表示

ちなみに、サムネイルはAdobe Fireflyで作ってみました。