FIBA W杯が開催されて日本代表の活躍が話題になりましたね。
私も現地までは行けませんでしたが、テレビの前で応援していました。
また強化試合は数試合に参加しました。 特にスロベニア戦は日本ベンチ裏の席がとれでテンションあがってましたね。
さて、W杯において日本代表がいたグループは、死のグループと言われていました。
確かに優勝したドイツ、NBA選手の方が多かったオーストラリア、勝てたとは言え長身揃いのフィンランド。
野暮ですが「日本代表って組み合わせ次第で決勝ラウンド行ける可能性あったんじゃね？」と思うわけです。
そこで、今回はデータを使って日本代表が決勝ラウンドに進出できたかどうかを簡単に検証していきたいと思います。
まあ、ネタ分析だと思ってください。大会側はくじ引きによる影響をなんとかしろとかいう意図はありません。
続きを読むtarget trial emulationについてまとめた中で不死時間バイアス:Immortal time biasを知りました。
Rのコードは以下の記事をかなりパクる参考にさせてもらっている。
フォローアップ開始から治療開始までに時間がある場合、その期間はイベントが起きない。 もし、イベントのが治療開始までにイベントが発生してしまった場合は、対照群に割り当てられることになる。 この時、フォローまでは絶対にイベントが起きないため、その分生存時間が伸びることになる。 このバイアスを不死時間バイアス：Immortal time biasと呼ぶ。
RでImmortal time biasが本当に発生するのかを確認する。
生存時間情報はワイブル分布を用いて生成する。
処置群のうち確率的に生成されたタイミングで処置が行われるとする。
この時、処置タイミングより前に生存時間に達してしまっている場合は、対象群に割り当てる。
一回1000サンプルを生成し、1000回繰り返すことで、バイアスの状況を分布として確認する。
library(tidyverse) library(survival) library(survminer) library("progress") gamma = 1.0 # ハザード比の設定 N = 1000 L = 1000 result_df <- data.frame() pb <- progress_bar$new(total = L) for(l in 1:L){ pb$tick() suvival_df <- data.frame(id = 1:N, Lambda = 1) %>% mutate(time1 = rweibull(N, shape=1.0, scale=Lambda), time2 = rweibull(N, shape=1.0, scale=Lambda*exp(-log(gamma))), treat = rbernoulli(N, p=0.5), treat_time = rweibull(N, shape=1.0, scale=0.6*Lambda), cons_time = runif(N, 0, 5) ) %>% mutate(time = if_else(treat==1, time2, time1), treat = if_else((time2 < treat_time)&(treat==T), F, treat), censoring = if_else((time > cons_time), 0, 1), time = pmin(time, cons_time)) fit_model <- coxph(Surv(time, censoring) ~ treat, data = suvival_df) result_df <- bind_rows(result_df, data_frame("i" = l, "Naive" = fit_model$coefficients[1])) }
結果を確認してみます。
result_longer_df <- result_df %>% pivot_longer(cols = -i, names_to = "col", values_to = "val") %>% mutate(HR = exp(val)) result_longer_df %>% ggplot() + aes(x=HR) + geom_histogram(colour="white", bins = 50) + geom_vline(aes(xintercept = gamma), colour="red") + theme_light() + xlab("HR") + ggtitle(label = str_c("Histograms of Estimates : True HR = ", round(gamma, 2)))
HR＝1と設定したのに、推定された処置群割り当ての影響はそれを下回る=イベント発生が過剰に抑制されている状態にあります。
> result_longer_df$HR %>% summary() Min. 1st Qu. Median Mean 3rd Qu. Max. 0.3892 0.5059 0.5450 0.5460 0.5819 0.8010
HR＝1.5と設定した場合でも同様でした。
> result_longer_df$HR %>% summary() Min. 1st Qu. Median Mean 3rd Qu. Max. 0.5077 0.6826 0.7290 0.7306 0.7766 1.0048
Immortal time biasの対応方法は複数提案されている。
処置群において処置が行われるまでの時間を削除して、処置後のみを対象とする方法。
ある時刻(landmerk)までの情報を削除して、処置後のみを対象とする方法。
最もImmortal time biasがない推定を行うことができるとされている。 時間依存性共変量を考えた生存時間分析を行う方法。TD-Cox modelとも呼ばれる。
具体的には処置が行われるまでが、一つのサンプルだと考える。つまり処置タイミングまでで打ち切りが行われるサンプルとする。 処置が行われている期間は、左側打ち切りだと考えて2つ目のサンプルとする。
Exclusion Method、Landmerk Method、Mantel-Byar Methodでどれほどバイアスが軽減されるかを確認する。
Landmerk Methodにおけるランドマークはt={0.3, 0.5}とした。
gamma_1 = 1.0 base_lambda = 0.6 N = 1000 L = 1000 result_df <- data.frame() pb <- progress_bar$new(total = L) for(l in 1:L){ pb$tick() suvival_df <- data.frame(id = 1:N, Lambda = base_lambda) %>% mutate(time1 = rweibull(N, shape=1.0, scale=Lambda), time2 = rweibull(N, shape=1.0, scale=Lambda*exp(-log(gamma_1))), treat = rbernoulli(N, p=0.7), treat_time = rweibull(N, shape=1.0, scale=0.6*Lambda), cons_time = runif(N, 0, 5) ) %>% rowwise() %>% mutate(time = if_else(treat==1, time2, time1), base_treat = if_else((time2 < treat_time)&(treat==T), F, treat), censoring = if_else((time > cons_time), 0, 1), time = pmin(time, cons_time)) # Mantel-Byar Method suvival_df2 <- tmerge(suvival_df, suvival_df, id = id, censoring = event(time, censoring)) suvival_df2 <- tmerge(suvival_df2, suvival_df2, id = id, treat = tdc(treat_time, base_treat)) %>% mutate(tstart = if_else(treat == FALSE & !is.na(treat), 0.0, tstart)) %>% group_by(id) %>% mutate(group = n()) %>% filter(!(is.na(treat) & base_treat == FALSE & group==2)) %>% replace_na(list("treat" = 0)) # Exclusion Method suvival_df3 <- suvival_df2 %>% filter(!((time > treat_time) & treat==FALSE)) %>% mutate(tstop = if_else((time > treat_time), tstop - treat_time, tstop), surv = Surv(tstop, censoring)) # landmerk method LM_time <- 0.5 suvival_df4_tmp <- suvival_df %>% filter(time >= LM_time) %>% mutate(time_event_LM1y = time - LM_time, time_treat_LM1y = treat_time - LM_time, flg_treat_LM_1yr = ifelse((treat_time >= LM_time | is.na(treat_time)), 0, 1), treat_LM1y = ifelse(time_treat_LM1y==1, 0,1) ) suvival_df4 <- tmerge(suvival_df4_tmp, suvival_df4_tmp, id = id, censoring = event(time, censoring)) suvival_df4 <- tmerge(suvival_df4, suvival_df4, id = id, treat = tdc(treat_time, treat_LM1y)) %>% replace_na(list("treat" = 0)) LM_time <- 1 suvival_df5_tmp <- suvival_df %>% filter(time >= LM_time) %>% mutate(time_event_LM1y = time - LM_time, time_treat_LM1y = treat_time - LM_time, flg_treat_LM_1yr = ifelse((treat_time >= LM_time | is.na(treat_time)), 0, 1), treat_LM1y = ifelse(time_treat_LM1y==1, 0,1) ) suvival_df5 <- tmerge(suvival_df5_tmp, suvival_df5_tmp, id = id, censoring = event(time, censoring)) suvival_df5 <- tmerge(suvival_df5, suvival_df5, id = id, treat = tdc(treat_time, treat_LM1y)) %>% replace_na(list("treat" = 0)) fit_nv <- coxph(Surv(time, censoring) ~ base_treat, data = suvival_df) fit_mb <- coxph(Surv(tstart,tstop, censoring) ~ treat, data = suvival_df2) #tstart,tstopと両方とも与える必要がある fit_em <- coxph(Surv(tstop, censoring) ~ treat, data = suvival_df3) fit_lm5 <- coxph(Surv(tstop, censoring) ~ treat, data = suvival_df4) fit_lm10 <- coxph(Surv(tstop, censoring) ~ treat, data = suvival_df5) result_df <- bind_rows(result_df, data_frame("i" = l, "Naive" = fit_nv$coefficients[1], "Mantel_Byar" = fit_mb$coefficients[1], "Exclusion" = fit_em$coefficients[1], "Landmerk5" = fit_lm5$coefficients[1], "Landmerk10" = fit_lm10$coefficients[1])) }
生存時間曲線を確認してみる。
result_longer_df <- result_df %>% pivot_longer(cols = -i, names_to = "col", values_to = "val") %>% mutate(val = exp(val)) ggplot(data=result_longer_df, aes(x=val)) + geom_histogram(colour="white", bins = 50) + facet_grid(col~.) + geom_vline(aes(xintercept = gamma_1), colour="red") + theme_bw() + ggtitle(label = str_c("Histograms of Estimates : HR = ", round(gamma_1, 2)))
Landmark MethodでもHR=1.0付近に分布しバイアスが小さく、Mantel Byar Methodが最も分散が小さい状態であった。 Landmark Methodはどれぐらいを削除するかを設定しなければならないことも考えると、Mantel Byar Methodを用いるのが良いと言える。
HR=0.7とHR=1.5の場合でも同様であった。
Target trial emulationという考え方を少し前に知りました。
面白そうなので、メモとしてまとめておきます。
(上記を中心とした参考資料の方が適切でわかりやすいのそちらを参照されたし)
Target trial emulationとは、観測データを用いてRCTの状態を再現することで適切な処置効果の推定を行おうというフレームワーク。
具体的な分析方法ではなくフレームワークであり、データを見る際により最適な状態に導いてくれる考え方がまとめられているものであるようだ。
Target trial emulationの考え方では、対象者の選択(Eligibility)、治療割付(Treatment Assignment)、タイムゼロの3つのタイミングが重要な要素となる。
最終的な分析対象を選ぶための条件が全て満たされるタイミング
処置を行ってたとしても、条件となる状況や行動等が実施されていない場合は、分析対象にはならない。
観測データでは、処置対象の状態(共変量)により、このタイミングが変化してしまうとバイアスが生じてしまうことに気をつけなければならない。
処理群に対して処置を行った(行い始めた)タイミング
観測データでは、処理割当が様々な共変量によりタイミングが変化することにある。
また、二重盲検試験とは異なり処置を担当する人物が、対象者の割り当て状況を理解しているため、間接的なバイアスが存在する可能性がある。
タイムゼロとは、フォローアップが開始されるタイミング
観測データの自然な開始は治療戦略が割り当てられたときであり、多くの場合、治療割当(処置割当)の開始の直前または開始と同時。
ランダム化後に開始すると、ランダム化とタイムゼロの間のすべてのサンプルが分析から除外されるため、選択バイアスが生じる可能性がある。
観察データを使用して、対象となる治験のタイムゼロをエミュレートする最も良い方法は、個人の治療が開始される時刻をタイムゼロと定義することです。
RCTの場合では「処置割付」と「タイムゼロ」が同じタイミングであり、「対象者選択」はそれらと同時かかすでに完了している状態と言える。 そして、これらの終了後に評価が行われている。
Target trial emulationでは、タイミングの順序がRCTと同様になるようなターゲットの選択や処理、分析を行うことになる。
タイムゼロから処置割当に時間があくと、処置割り当てまでに対象者がフォローアップから外れる(何かしらの理由で観測が中止される/死亡する)ことになる。 処置割当までフォローアップできたということにより、生存時間が比較的長くなったり、差が生じてしまうことがある。 このようなバイアスをimmortal time bias(不死時間バイアス)と呼ばれている。
フォローアップ開始から処置開始までに時間がある場合、その期間はイベントが起きない。 もし、イベントが処置開始までに発生してしまった場合は、対照群に割り当てられることになる。 この時、処理群ではフォローまでは絶対にイベントが起きない状態になるため、その分生存時間が伸びてしまう場合がある。
いくつかの対処方法が考えられるが、生存時間分析の場合は処置の有無を時間依存変量として扱う方法があげられる。 これについて、次のブログ記事で扱いたい。
処置割付のタイミングが曖昧であったり、データの取得方法によってずれてしまう場合がある。
例えば、特定のイベントを処置だと考えると、そのイベントが直後からのフォローアップが難しい場合などが考えられる。
この時にフォローアップが実施される割合が処置によって変化したり、フォローアップの実施と未実施で結果が変化してしまう場合に、選択バイアスが生じてしまう。
また、フォローアップを開始しようとした時にすでに処置を行っているような対象が存在する場合、処置までにフォローアップが終了等により処置が実施されないことになる。
このような場合に発生するバイアスのことをPrevalent user biasと呼ばれる。
例えば、ある薬の効果を確認する時すでにその薬を服用している人がいる可能性があり、その人物を処置群に割り当てる場合などが挙げられる。 このバイアスもimmortal time biasと同様に、処置を開始するときにその薬を服用しているということは、その時点までその薬の服用をやめるようなイベント(治療の停止/死亡)などが起こっていないために発生する。 処置群には生き残っている人たちが多く含まれることになり、対照群よりも生存時間が過剰に伸びてしまうことになる。
末期腎不全患者に対する維持透析開始時の腎機能と予後の関係を検討したスウェーデンの観察研究
観測データを通常の解析を行うと、透析開始時点と腎機能の予後に有意な関係性は見られなかった。
しかしこの研究では、以下の二点が問題であった。
タイムゼロは透析開始時ではなくある条件と満たした時点であり、選択バイアスが発生してしまう。比較的健康的な患者は透析導入が遅くなる傾向がある。
治療割当は観測開始時点とした場合、透析開始まだ亡くならなかった患者が対象にあり、不死時間バイアスが発生していると考えられる。遅く透析を導入した患者は、長生きするような傾向に見えてしまう。
それぞれの対応方法として、タイムゼロは条件と満たしたと時点とし、治療割当は過去の治療計画情報を含めた期間も解析対象にしている。
結果として、4年早く透析を導入すると1.6カ月の生存期間が伸びることが示された。
参考
immortal time biasが存在する場合の生存時間分析はいくつかの資料で取り上げられていたが、それ以外のパターンの対応方法で具体的なもの/汎用的なものは少なかった。 対応方法は個別に検討していかなければならず、実際の勘所がまだよくわからない。 (Target trial emulationはフレームワークに過ぎないので、仕方ないとは思うが...)
