【学習動機】
Immortal time bias (不死時間バイアス)というバイアスがある. フォローアップ開始から治療開始までに時間がある場合, その分の時間はイベントが起きないので (もしもイベントが起きていたら, 未治療群のイベントとして扱われる), その分の時間(イベントが死亡であるならば「不死」の時間となる)が起因となってバイアスが 生じる. immortal time biasへの対処法を学んでいたら, Mi et al. 2016 “Use of the landmark method to address immortal person-time bias in comparative effectiveness research: a simulation study”を見つけた. シミュレーション研究でimmortal time biasを含むデータに対して, naiveな方法, 未治療時間を削除する方法, ランドマーク解析, Mantel-Byarの方法を比較している. 論文ではSASが使われている(コードはない)が今回はこれをRで再現してみる.
【学習内容】
基本設定とパラメータ設定
パッケージ読み込み.
# パッケージ読み込み ########################
library(tidyverse)
library(survival)
library(survminer)サンプルサイズは1000として(論文では10000だが速度の都合で1000としている), フォローアップ期間は5年間(5*365日, うるう年は簡単化のため無視)とする.
その他, 必要なパラメータ設定. 未治療群に対する治療群のハザード比(HR)は0.7と設定する.
# 基本設計 #################################
set.seed(2022)
N <- 1000 # sample size(論文では10000だが速度の都合で1000としている)
fup_max <- 5 * 365 # 5years, max followup
# パラメータ設定 #############################
alpha_0 <- -1.14
alpha_L <- log(1.1)
alpha_M <- log(1.25)
alpha_H <- log(1.5)
# Hazard Ratio(HR)
HR_0 <- 0.7 # 0.5, 0.7, 1.0, 1.3, 1.5の変更オプション
gamma_0 <- log(HR_0) # log HR
gamma_0## [1] -0.3566749# beta
# beta < 1 で時間経過で治療確率は減少
# beta = 1 で時間経過で治療確率は一定
# beta > 1 で時間経過で治療確率は増加
# treatment timeのパラメータ
beta_T <- 1 # 0.5, 1, 1.5の変更オプション
lambda_T <- 1/800 # 1/800
# event timeのパラメータ
beta_E <- 1 # 0.5, 1, 1.5の変更オプション
lambda_E <- 1/500データ生成
設定したパラメータを使って, シミュレーションデータを生成していく.
(修正(2022-07-01) : twitterでnissinboさんからコメントをいただき,
treatの変数の書く位置を間違えていたことに気付けたので訂正しています. コメント, ありがとうございました. Naive MethodとExclusion Methodの結果が修正で変わっています.)
# データ生成 #############################
df0 <- data.frame(
# 共変量 X1~X8 の生成
X1 = sample(x = c(1,0), size = N, replace = TRUE, prob = c(0.3,0.7)),
X2 = sample(x = c(1,0), size = N, replace = TRUE, prob = c(0.4,0.6)),
X3 = sample(x = c(1,0), size = N, replace = TRUE, prob = c(0.5,0.5)),
X4 = sample(x = c(1,0), size = N, replace = TRUE, prob = c(0.5,0.5)),
X5 = rnorm(N, mean=0, sd=1),
X6 = rnorm(N, mean=0, sd=1),
X7 = rnorm(N, mean=0, sd=1),
X8 = rnorm(N, mean=0, sd=1))
df1 <- df0 %>%
mutate(
c0 = exp(alpha_0 + alpha_L*X1 + alpha_M*X2 + alpha_H*X3 + alpha_L*X5 + alpha_M*X6 + alpha_H*X7),
c1 = exp(alpha_0 + alpha_L*X1 + alpha_M*X2 + alpha_H*X4 + alpha_L*X5 + alpha_M*X6 + alpha_H*X8),
# 逆変換に使用する一様乱数の生成
u0 = runif(N, min=0, max=1),
u1 = runif(N, min=0, max=1)
)
# treatment(exposure) timeの乱数生成 : ワイブル分布
ft0 <- function(u0, lambda0, beta0, c0){
1/lambda0 * ((-log(1-u0))/c0)^(1/beta0)
}
# event timeの乱数生成(場合分け) : 時間依存性共変量の治療Zに対応
ft1_1 <- function(u1, lambda1, beta1, c1){
1/lambda1 * ((-log(1-u1))/c1)^(1/beta1)
}
ft1_2 <- function(u2, lambda2, beta2, b2, c1, c2){
((-log(1-u2)-c1*(lambda2*b2)^(beta2))/(c2*lambda2^beta2) + b2^beta2)^(1/beta2)
}
df2 <- df1 %>%
mutate(
t0 = ft0(u0=u0, lambda0 = lambda_T, beta0 = beta_T, c0 = c0), # treatment time
Z = ifelse(t0<=fup_max, 1, 0), # 1:treated, 0:untreated
#treat = ifelse(Z==1, "yes", "no"), (ミス. eventの時間との兼ね合いがあるので, ここで治療はまだ定義できない)
c2 = exp(alpha_0 + alpha_L*X1 + alpha_M*X2 + alpha_H*X4 + alpha_L*X5 + alpha_M*X6 + alpha_H*X8 + gamma_0*Z),
t1_1 = ft1_1(u1 = u1, lambda1 = lambda_E, beta1 = beta_E, c1 = c1),
t1_2 = ft1_2(u2 = u1, lambda2 = lambda_E, beta2 = beta_E, b2 = t0, c1 = c1, c2 = c2)
)
df3 <- df2 %>%
mutate(
t1 = ifelse(u1<1-exp(-c1*(lambda_E*t0)^beta_E), t1_1, t1_2), # event time
time_event = ifelse(t1 < fup_max, t1, fup_max), # event time or max followup
time_treat = ifelse(t0 < fup_max & t0 < time_event, t0, NA), # treat tiem
flg_treat = ifelse(t0 < fup_max & t0 < time_event, 1, 0), # 1:treated, 0:untreated
treat = ifelse(flg_treat==1, "yes", "no"), # こっちに修正
flg_event = ifelse(t1 < fup_max, 1, 0)
)推定方法
methodの前のABCDは論文に対応. 実装しやすい順番でコードを書いている. ここでは C : Naive Method, B : Exclusion Method, D : Landmark Method, A : Mantel-Byar Method の順で書いている.
カプランマイヤー曲線を書いて, Coxモデルを当てはめている. 共変量は正しくモデル化できているという設定.
C : Naive Method
何もせず, ナイーブに解析する方法.治療群の方が長生きするはずだが, ほとんど変わらない結果になっている.
予想通り, 治療群の方が随分長生きする結果になっている.
# C : Naive Method############################################
df_naive <- df3
res_KM_naive <- survival::survfit(Surv(time_event, flg_event)~treat, data=df_naive)
survminer::ggsurvplot(res_KM_naive, dat=df_naive, legend.title = "legend",
legend=c(.8,.8), title = "KM curve ( Naive Method )",
legend.labs=c('untreated','treated'), size=1,
xlab="days", ylab="Event free survival Rate",
censor=T, ggtheme = theme_gray())
res_cox_naive <- survival::coxph(Surv(time_event, flg_event)~treat+X1+X2+X4+X5+X6+X8, data = df_naive)
res_cox_naive## Call:
## survival::coxph(formula = Surv(time_event, flg_event) ~ treat +
## X1 + X2 + X4 + X5 + X6 + X8, data = df_naive)
##
## coef exp(coef) se(coef) z p
## treatyes -1.26784 0.28144 0.09149 -13.858 < 2e-16
## X1 0.13547 1.14507 0.08031 1.687 0.091642
## X2 0.26848 1.30798 0.07575 3.544 0.000394
## X4 0.45492 1.57606 0.07584 5.998 1.99e-09
## X5 0.06121 1.06312 0.03806 1.608 0.107750
## X6 0.32148 1.37916 0.03807 8.445 < 2e-16
## X8 0.40087 1.49312 0.03938 10.180 < 2e-16
##
## Likelihood ratio test=434.6 on 7 df, p=< 2.2e-16
## n= 1000, number of events= 727B : Exclusion Method
治療された場合, 未治療の期間を単純に除外する方法.治療群の方が長生きするはずだが, むしろ未治療の方が長生きする結果になっている.
予想通り, 治療群の方がやや長生きする結果になっている.
# B : Exclusion Method############################################
df_EM <- df3 %>%
mutate(
time_event_EM = ifelse(flg_treat==1, time_event-time_treat, time_event)
)
res_KM_EM <- survival::survfit(Surv(time_event_EM, flg_event)~treat, data=df_EM)
survminer::ggsurvplot(res_KM_EM, dat=df_EM, legend.title = "legend",
legend=c(.8,.8), title = "KM curve ( Exclusion Method )",
legend.labs=c('untreated','treated'), size=1,
xlab="days", ylab="Event free survival Rate",
censor=T, ggtheme = theme_gray())
res_cox_EM <- survival::coxph(Surv(time_event_EM, flg_event)~treat+X1+X2+X4+X5+X6+X8, data = df_EM)
res_cox_EM## Call:
## survival::coxph(formula = Surv(time_event_EM, flg_event) ~ treat +
## X1 + X2 + X4 + X5 + X6 + X8, data = df_EM)
##
## coef exp(coef) se(coef) z p
## treatyes -0.78773 0.45488 0.09049 -8.705 < 2e-16
## X1 0.14052 1.15088 0.08042 1.747 0.080576
## X2 0.27205 1.31265 0.07574 3.592 0.000328
## X4 0.41540 1.51498 0.07567 5.490 4.03e-08
## X5 0.05902 1.06079 0.03809 1.549 0.121283
## X6 0.29858 1.34794 0.03807 7.843 4.40e-15
## X8 0.38027 1.46269 0.03960 9.603 < 2e-16
##
## Likelihood ratio test=279.5 on 7 df, p=< 2.2e-16
## n= 1000, number of events= 727D : Landmark Method
ランドマークを1年に設定したものと2年に設定したものを作成.
どちらの結果も, 想定通り, 治療群の方が良い結果になっている.
# D : Landmark Method############################################
LM_year <- 1
df_LM_1yr <- df3 %>%
filter(time_event >= 365*LM_year) %>%
mutate(
time_event_LM_1yr = time_event - 365*LM_year,
time_treat_LM_1yr = time_treat - 365*LM_year,
flg_treat_LM_1yr = ifelse((time_treat >= 365*LM_year | is.na(time_treat)), 0, 1),
treat_LM_1yr = ifelse(flg_treat_LM_1yr==1, "yes","no")
)
res_KM_LM_1yr <- survival::survfit(Surv(time_event_LM_1yr, flg_event)~treat_LM_1yr, data=df_LM_1yr)
survminer::ggsurvplot(res_KM_LM_1yr, dat=df_LM_1yr, legend.title = "legend",
legend=c(.9,.9), title = "KM curve ( Landmark Method (1yr) )",
legend.labs=c('untreated','treated'), size=1,
xlab="days", ylab="Event free survival Rate",
censor=T, ggtheme = theme_gray())
res_cox_LM_1yr <- survival::coxph(Surv(time_event_LM_1yr, flg_event)~treat_LM_1yr+X1+X2+X4+X5+X6+X8, data = df_LM_1yr)
res_cox_LM_1yr## Call:
## survival::coxph(formula = Surv(time_event_LM_1yr, flg_event) ~
## treat_LM_1yr + X1 + X2 + X4 + X5 + X6 + X8, data = df_LM_1yr)
##
## coef exp(coef) se(coef) z p
## treat_LM_1yryes -0.35527 0.70098 0.12564 -2.828 0.00469
## X1 0.15553 1.16828 0.10372 1.499 0.13375
## X2 0.12633 1.13466 0.09690 1.304 0.19232
## X4 0.56087 1.75219 0.09597 5.844 5.09e-09
## X5 0.04396 1.04495 0.04841 0.908 0.36375
## X6 0.24623 1.27920 0.04864 5.062 4.14e-07
## X8 0.41710 1.51755 0.04986 8.365 < 2e-16
##
## Likelihood ratio test=123.6 on 7 df, p=< 2.2e-16
## n= 724, number of events= 451LM_year <- 2
df_LM_2yr <- df3 %>%
filter(time_event >= 365*LM_year) %>%
mutate(
time_event_LM_2yr = time_event - 365*LM_year,
time_treat_LM_2yr = time_treat - 365*LM_year,
flg_treat_LM_2yr = ifelse((time_treat >= 365*LM_year | is.na(time_treat)), 0, 1),
treat_LM_2yr = ifelse(flg_treat_LM_2yr==1, "yes","no")
)
res_KM_LM_2yr <- survival::survfit(Surv(time_event_LM_2yr, flg_event)~treat_LM_2yr, data=df_LM_2yr)
survminer::ggsurvplot(res_KM_LM_2yr, dat=df_LM_2yr, legend.title = "legend",
legend=c(.9,.9), title = "KM curve ( Landmark Method (2yr) )",
legend.labs=c('untreated','treated'), size=1,
xlab="days", ylab="Event free survival Rate",
censor=T, ggtheme = theme_gray())
res_cox_LM_2yr <- survival::coxph(Surv(time_event_LM_2yr, flg_event)~treat_LM_2yr+X1+X2+X4+X5+X6+X8, data = df_LM_2yr)
res_cox_LM_2yr## Call:
## survival::coxph(formula = Surv(time_event_LM_2yr, flg_event) ~
## treat_LM_2yr + X1 + X2 + X4 + X5 + X6 + X8, data = df_LM_2yr)
##
## coef exp(coef) se(coef) z p
## treat_LM_2yryes -0.32646 0.72148 0.13701 -2.383 0.017187
## X1 0.24612 1.27905 0.13767 1.788 0.073815
## X2 0.22257 1.24928 0.12858 1.731 0.083451
## X4 0.59280 1.80905 0.12836 4.618 3.87e-06
## X5 -0.03450 0.96609 0.06531 -0.528 0.597326
## X6 0.21789 1.24345 0.06434 3.387 0.000708
## X8 0.37763 1.45882 0.06738 5.604 2.09e-08
##
## Likelihood ratio test=65.21 on 7 df, p=1.367e-11
## n= 527, number of events= 254A : Mantel-Byar Method
時間依存性共変量を考慮したモデルという理解で実装(Mantel-Byar methodがこれを指しているという理解だが, 自信はない ) こちらも想定通り, 治療群の方が良い結果になっている.
# A : Mantel-Byar Method############################################
df3_MB1 <- df3 %>% mutate(id=1:nrow(df3))
df3_MB2 <-
df3_MB1 %>% tmerge(df3_MB1, id=id, flg_event_MB=event(time_event,flg_event)) %>%
tmerge(df3_MB1, id=id, flg_treat_MB=tdc(time_treat))
res_KM_MB <- survival::survfit(Surv(tstart,tstop,flg_event_MB)~flg_treat_MB, data=df3_MB2)
#plot(res_KM_MB, col=c("red","blue"))
#legend("topright", names(res_KM_MB$strata),lty=1,col=c("red","blue"))
survminer::ggsurvplot(res_KM_MB, dat=df3_MB2, legend.title = "legend",
legend=c(.9,.9), title = "KM curve ( Mantel-Byar Method )",
legend.labs=c('untreated','treated'), size=1,
xlab="days", ylab="Event free survival Rate",
censor=T, ggtheme = theme_gray())
res_cox_MB <- survival::coxph(Surv(tstart,tstop,flg_event_MB)~flg_treat_MB+X1+X2+X4+X5+X6+X8, data = df3_MB2)
res_cox_MB## Call:
## survival::coxph(formula = Surv(tstart, tstop, flg_event_MB) ~
## flg_treat_MB + X1 + X2 + X4 + X5 + X6 + X8, data = df3_MB2)
##
## coef exp(coef) se(coef) z p
## flg_treat_MB -0.45680 0.63330 0.09665 -4.727 2.28e-06
## X1 0.19576 1.21623 0.08013 2.443 0.01457
## X2 0.20292 1.22497 0.07565 2.682 0.00731
## X4 0.52038 1.68267 0.07575 6.870 6.44e-12
## X5 0.04905 1.05027 0.03803 1.290 0.19715
## X6 0.26768 1.30692 0.03826 6.997 2.62e-12
## X8 0.44220 1.55612 0.03921 11.277 < 2e-16
##
## Likelihood ratio test=232.5 on 7 df, p=< 2.2e-16
## n= 1343, number of events= 727推定結果の比較
真の値と各推定値を比較してみる. ランドマーク解析法とMantel-Byarの方法が近いところを推定している結果になった.
とくに, 1年のランドマーク解析はかなり近いところを推定している.
# 推定結果の比較(Log HR)
data.frame(
LHR_true = gamma_0,
LHR_naive = res_cox_naive$coefficients[1],
LHR_EM = res_cox_EM$coefficients[1],
LHR_LM_1yr = res_cox_LM_1yr$coefficients[1],
LHR_LM_2yr = res_cox_LM_2yr$coefficients[1],
LHR_MB = res_cox_MB$coefficients[1]
) %>% t() %>% as.data.frame() %>%
rename(val=treatyes) %>%
mutate(exp=exp(val))## val exp
## LHR_true -0.3566749 0.7000000
## LHR_naive -1.2678361 0.2814400
## LHR_EM -0.7877292 0.4548765
## LHR_LM_1yr -0.3552714 0.7009832
## LHR_LM_2yr -0.3264568 0.7214756
## LHR_MB -0.4568048 0.6333040【学習予定】
シミュレーション1回分を実装できた(と思う)ので, 次回でこのシミュレーションをたくさん繰り返すコードを書いて実験する. その2に続く.
