Immortal time bias

不死の時間とは、デザイン上、死亡やイベントが発生しない追跡期間のこと。

例えば、退院日が追跡調査の開始日となる場合に、退院後に処方箋が調剤されるのを待つといったことです。
この待ち時間が存在するということはその期間は死んでいないことを意味し、治療群または被曝群となる個人は、治療の定義が満たされるまで生存していなければならないからです。治療を受ける前にイベントが発生した場合は、未治療または未曝露のグループに属することになる。この「不死の期間」が治療状況に関して誤って分類されたり、解析から除外されたりすると、バイアスが生じる。

不死の期間バイアスは、治療群に偽の生存率の優位性を与えることにより、必然的に研究対象の治療に有利な結果に偏るため、特に問題となる。

対処

対処は3つ
- 時点マッチング - ランドマーク解析
- 時間依存性解析

今回は時間依存性解析をRで行う。

Time varying covariate

時間依存供変量とは時間とともにその値を変化させるような変数
2つに分類される

内的時間依存性供変量:
投薬開始や心臓移植などの個人内要因による変化
- スタンフォード心臓移植
- アルコール肝硬変のPTの経時変化と生存率

外的時間依存供変量:
- 大気汚染など外的要因で変化
- スウェーデンの小麦の値段と５０歳以上の生存率

Time variyng covariate in Cox model

Coxモデルの長所の1つは、時間とともに変化する共変量を包含できること

時間依存の共変量が機能する理由は、Coxモデルの基本的な動作方法に基づく

Coxでは人ではなく人ー時間で、イベントの発生確率を考える。

各イベント発生時に、 イベントが発生した被験者の現在の共変量の値と
その時点でリスクがあったすべての被験者の現在の値を比較する

library

では実際のデータを使ってやってみる。 データはスタンフォードの心臓移植のデータを使用する。 研究仮説は

“移植により予後が改善する”

である。

このデータはsurvivalの中に入っている。 データの概要は?jasaで確認のこと。

library(survival)
library(tidyverse)
library(rms)
library(broom)

jasaデータから必要な変数を選択する。 先に、immortal timeに対応しない場合にどうなるのかを見ておく。

dat<-jasa %>% select(accept.dt,tx.date,fu.date,fustat,transplant,surgery)
dat$id <- 1:nrow(dat) # idを付与

dat$futime <- pmax(0.5,dat$fu.date-dat$accept.dt) #組入から終了までのfollow up timeを作成
dat$txtime <- ifelse(dat$tx.date==dat$fu.date,
                     dat$futime-0.5,dat$tx.date-dat$accept.dt) #暴露までの期間を作成

#データの都合上、0になるセルがでるので0.5を足しておく。
  
head(dat)

futimeが生存期間、fustatが転帰上の状態で1が死亡、0が生存なので、これらから生存時間データを作成し、カプランマイヤー曲線を書く

surv.tv<-with(dat,Surv(futime,fustat))
surv.km<-survfit(surv.tv~transplant,data=dat)
#plot(surv.km, col=(c("black","red"))) #plotでもかけるけど今回は下で

library(survminer)
##  要求されたパッケージ ggpubr をロード中です
ggsurvplot(surv.km,dat=dat,size=1, pval=T,size=1,conf.int = T,
           risk.table = T, palette="Lancet")

ハザード比を求めてみる。

f_bias<-coxph(surv.tv~transplant,data=dat)
tidy(f_bias, exponentiate = T, conf.int = T)

HR0.27 (95%CI 0.17-0.43)という強烈な効果を示している

はたしてこれはガチの効果か、バイアスか。

bias

ですね。

immortal time biasを可視化してみます。

dat %>% group_by(transplant) %>% 
  ggplot(aes(x=reorder(reorder(id,futime),transplant),y=futime, 
             fill=transplant, group=transplant))+
  geom_col(stat = "identity", position = "dodge")+
  coord_flip()
## Warning: Ignoring unknown parameters: stat

青が移植群の生存時間
黒が移植なし群の生存時間です。
圧倒的に、移植有りのほうが長そうですが、待機時間が含まれています。

それを可視化してみます。

dat %>% group_by(transplant) %>% 
  ggplot(aes(x=reorder(reorder(id,futime),transplant),fill=transplant,
             group=transplant))+
  geom_col(aes(y=futime),stat = "identity", position = "dodge")+
  geom_col(aes(y=txtime,col="waiting time"),stat = "identity", position = "dodge")+
  coord_flip()+
  theme_bw()+ #枠線と目盛り線のあるテーマを設定
  theme(panel.grid = element_blank()) #目盛り線を消去
## Warning: Ignoring unknown parameters: stat

## Warning: Ignoring unknown parameters: stat
## Warning: Removed 34 rows containing missing values (geom_col).

赤い時間がImmortal time biasです。 この赤い期間が、率比の分母に含まれるので治療群で過剰に率比が小さくなります。

\(Rate_{treated}=(治療群の死亡)/(死亡までの治療期間の合計人-時間)\)
\(Rate_{untreated}=(非治療群の死亡)/(死亡までの非治療期間の合計人-時間)\)

\(Rate Ratio=Rate_{treated}/Rate_{untreated}\)

これがimmortal time biasです。

移植の待機時間のせいで治療である移植が時間依存性に0→1に変化するために起こっています。

time dependent covariate

共変量「移植」を時間で分類するためには、実際には次のようなデータセットが必要です。

head(jasa1)
# head(heart) #heart data is maybe better?

このデータが実際に解析に用いられる縦型のデータです。

ID1と2は移植されず、停止時間に死亡。 ID3は初日に移植され、その後死亡。 ID4は35日目に移植され、38日目に死亡。 ID4は2列に分かれているので、1列目のデータは非移植、2列目のデータは移植されたものとして扱われます。これが時間依存共変量解析の簡単な例です。

では、jasa dataをjasa1に変換するためのコードを見てみましょう。

悲しいことに、コホートエントリーと同じ日に死亡した被験者がいて、開始-停止時間が(0,0)となっていますが、これは無効なので、死亡した日が0.5日になるようにデータを修正しなければなりません。

また、移植と同日に死亡した被験者もいますが、これは移植群での死亡とみなすべきなので、移植日を0.5日後ろにずらす必要があります。

head(dat)
dat<-jasa %>% select(accept.dt,tx.date,fu.date,fustat,transplant,surgery)
dat$id <- 1:nrow(dat)
dat$futime <- pmax(0.5,dat$fu.date-dat$accept.dt)
dat$txtime <- ifelse(dat$tx.date==dat$fu.date,
                     dat$futime-0.5,dat$tx.date-dat$accept.dt)

そして、survival::tmerge()を使います。

newdat <- tmerge(dat,dat,id=id,event=event(futime,fustat))
newdat[1:5, 7:12]

eventはtmergeのオプションです。 newname = event(y,x) 時刻yのイベントをマークします。
xが欠落している通常のケースでは、新しい0/1変数は、生存時間の0/1ステータス変数と同様になります。

newdat <- tmerge(newdat,dat,id=id,tx=tdc(txtime))

newdat[1:5, 7:13]

newname = tdc(y, x, init)により新しい時間依存の共変量変数が作成されます。

引数yは、開始時間と終了時間の尺度であると仮定され、各インスタンスはその時間における「条件」の発生を記述します。第2引数のxはオプションです。xがない場合、カウント変数は各被験者の0から始まり、イベントの発生時に1になります。xが存在する場合、時間依存性共変量の値は、存在すればinitの値、存在しなければtdcstartオプションの値に初期化され、各観測でxの値に更新されます。オプション na.rm=TRUE の場合、x の欠損値はまず除去されます。

これでjasa1のデータとおなじになりました

head(jasa1)

ではこのデータを元に、時間依存性Coxモデルを試してみます。

newdat$Surv <- Surv(newdat$tstart,newdat$tstop,newdat$event)
plot(survfit(Surv~tx,data=newdat),col=c("red","blue"))

ggsurvplot(survfit(Surv~tx,data=newdat), conf.int = T, risk.table = T, palette="Lancet")

ハザード比を求めます。 最初はクラスターなしで求めます。

f1<-coxph(Surv~tx+surgery,data=newdat)
summary(f1)
## Call:
## coxph(formula = Surv ~ tx + surgery, data = newdat)
## 
##   n= 170, number of events= 75 
## 
##            coef exp(coef) se(coef)      z Pr(>|z|)  
## tx       0.1562    1.1690   0.2965  0.527   0.5984  
## surgery -0.7491    0.4728   0.3596 -2.083   0.0373 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## tx         1.1690     0.8554    0.6538    2.0903
## surgery    0.4728     2.1150    0.2337    0.9567
## 
## Concordance= 0.563  (se = 0.031 )
## Likelihood ratio test= 5.36  on 2 df,   p=0.07
## Wald test            = 4.54  on 2 df,   p=0.1
## Score (logrank) test = 4.72  on 2 df,   p=0.09

ID共通者では相関があるはずなのでクラスターしてみる。

newdat
f<-coxph(Surv~tx+surgery+cluster(id),data=newdat)
summary(f)
## Call:
## coxph(formula = Surv ~ tx + surgery, data = newdat, cluster = id)
## 
##   n= 170, number of events= 75 
## 
##            coef exp(coef) se(coef) robust se      z Pr(>|z|)  
## tx       0.1562    1.1690   0.2965    0.2975  0.525    0.600  
## surgery -0.7491    0.4728   0.3596    0.3387 -2.211    0.027 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## tx         1.1690     0.8554    0.6525    2.0944
## surgery    0.4728     2.1150    0.2434    0.9184
## 
## Concordance= 0.563  (se = 0.032 )
## Likelihood ratio test= 5.36  on 2 df,   p=0.07
## Wald test            = 5  on 2 df,   p=0.08
## Score (logrank) test = 4.72  on 2 df,   p=0.09,   Robust = 4.61  p=0.1
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

クラスターにしても変わらなかった。

なお個人ごとに複数のイベントがある場合を除いて、個人が複数の行を持つ可能性があることを分析で考慮する必要はない。 尤度方程式では、個人の時間間隔は重ならないので、どの時点でも、個人ごとに最大で1つの行の情報しか使用しないため
今回は複数イベントでないのでクラスターしなくても結果は変わらないよう

weak point

実際に時間依存性Coxをやってみて思ったのは、本当にこれで良いのか感。

  • 観察前のImmortal timeには対応しない（selection biasを補正できない）

  • Time varying Cox modelではUntreated側の分母に組み込むため保守的になる？問い合わせ中

  • 対照群扱いになるImmortal time部分のデータは実際には早期に打ち切り扱いになっている。これはいいのか？選択バイアスになっていないか気になるところ。

これならUnbiasだよと教わるけど…。誰か分かる人は教えて下さい