IPTW 競合リスク累積発生率曲線の書き方

IPTW で重み付けした競合リスク累積発生曲線を書く方法

＞＞もう統計で悩むのを終わりにしませんか？

↑1万人以上の医療従事者が購読中

IPTW 競合リスク累積発生率曲線の書き方の前に IPTW、競合リスク、累積発生曲線とは？
IPTW 競合リスク累積発生率曲線の書き方
IPTW 競合リスク累積発生率曲線の書き方 IPTW の計算
IPTW 競合リスク累積発生率曲線の書き方累積発生率の推定
IPTW 競合リスク累積発生率曲線実際の描画
まとめ
参考 PDF

IPTW 競合リスク累積発生率曲線の書き方の前に IPTW、競合リスク、累積発生曲線とは？

IPTW に関しては以下の過去記事を参照

toukeier.hatenablog.com

競合リスクに関しては以下の過去記事を参照

toukeier.hatenablog.com

累積発生率（累積イベント）曲線に関しては以下の過去記事を参照

toukeier.hatenablog.com

IPTW 競合リスク累積発生率曲線の書き方

R の cmprskcoxmsm パッケージを使う方法を解説する

install.packages('cmprskcoxmsm')
library(cmprskcoxmsm)

サンプルデータは同じパッケージに含まれている follic というデータ

データセットには以下の変数が入っている

age 年齢
hgb ヘモグロビン値（ g/L ）
clinstg 臨床ステージ： 1= ステージ I 、2= ステージ II
ch 化学療法： Y= 実施あり、N= 実施なし
rt 放射線療法：全員 Y= 実施ありのため今回は使わず
time 最初のイベントが起きた時間
status イベントの理由： 1= 再燃、2= 死亡、0= 無効

＞＞もう統計で悩むのを終わりにしませんか？

↑1万人以上の医療従事者が購読中

IPTW 競合リスク累積発生率曲線の書き方 IPTW の計算

IPTW の計算には、doPS() 関数を使う

データ、治療（二値のみ。ここでは化学療法 ch）、治療のありの値、IPTW を作成するための交絡因子を引数として指定する

D <- doPS(follic, Trt='ch', Trt.name='Y', 
VARS.=c('age', 'hgb', 'clinstg'))

IPTW の計算結果が含まれる、D の先頭 5 行を表示させると以下のとおりである

ps_ate が傾向スコアで、ipw_ate_unstab が傾向スコアの逆数をとるだけのいわゆる IPTW、ipw_ate_stab が群の割合をかけた安定化 IPTW である

IPTW で交絡因子がどのくらい似通った分布になっているかを確認するために、plot.PS() 関数を実行する

plot.PS(D)

すると以下のグラフが描かれる

交絡因子と考えた、age, hgb, clinstg の３つとも、weighted のほうが Standard Mean Difference にして -0.1 から 0.1 の範囲に収まっていて、重みづけはうまくいっている

IPTW 競合リスク累積発生率曲線の書き方累積発生率の推定

次に、cif_est() 関数を使って、IPTW 累積発生率を推定する

データ、イベントまでの時間、イベントの変数名、イベントの値の種類、どのイベントに注目するか、重みづけの種類を指定する

CIF <- cif_est(data=D$Data, time='time', Event.var='status', 
  Events=c(0,1,2), 
  cif.event=1,
  weight.type='Unstabilized')

この推定には数分時間がかかった（CPU: Ryzen 3 3300X 4 コア 8 スレッド）

出来上がりのデータセットの先頭 5 行は以下の通りである

この中の、time, cif.control（化学療法なし）, cif.expsure（化学療法あり）を使って、IPTW 累積発生率曲線を書いていく

IPTW 競合リスク累積発生率曲線実際の描画

ggplot() 関数を使って書いていく

スクリプトは以下の通り

ggplot(CIF$cif_data, aes(x=time, y=cif.control, 
    color='No Chemotherapy'))+
  geom_line(linewidth=1.4)+
  geom_line(aes(x=time, y=cif.exposure, 
    color='Chemotherapy'), linewidth=1.4)+
  labs(x='Time', y='IPTW Cumulative incidence functuon')+
  ylim(c(0,1))