【学習動機】
その1のつづき.
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”のシミュレーション(naiveな方法, 未治療時間を削除する方法, ランドマーク解析(Landmark analysis), Mantel-Byarの方法を比較している.)をRでreplicateしてみている.
【学習内容】
パッケージ読み込み.
# パッケージ読み込み ########################
library(tidyverse)
library(survival)
library(survminer)前回のコードを活用して, シミュレーション回数と各シミュレーションのサンプルサイズと真のハザード比を決定して, シミュレーションを回し, 推定値のヒストグラムとバイアスとMSEが出力される関数を作成.
(修正(2022-07-01) : twitterでnissinboさんからコメントをいただき,
treatの変数の書く位置を間違えていたことに気付けたので訂正しています. コメント, ありがとうございました. Naive MethodとExclusion Methodの結果が修正で変わっています.)
# 関数化 ######################################################
sim_Mi2016 <- function(nsim, N, HR_0){
gamma_0 <- log(HR_0)
# 基本設計 #################################
fup_max <- 5 * 365 # 5years, max followup
LHR_true <- as.numeric(nsim)
LHR_naive <- as.numeric(nsim)
LHR_EM <- as.numeric(nsim)
LHR_LM_1yr <- as.numeric(nsim)
LHR_LM_2yr <- as.numeric(nsim)
LHR_MB <- as.numeric(nsim)
# パラメータ設定 #############################
alpha_0 <- -1.14
alpha_L <- log(1.1)
alpha_M <- log(1.25)
alpha_H <- log(1.5)
# 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
for(i in 1:nsim){
# データ生成 #############################
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)
)
# C : Naive Method############################################
df_naive <- df3
res_cox_naive <- survival::coxph(Surv(time_event, flg_event)~treat+X1+X2+X4+X5+X6+X8, data = df_naive)
# B : Exclusion Method############################################
df_EM <- df3 %>%
mutate(
time_event_EM = ifelse(flg_treat==1, time_event-time_treat, time_event)
)
res_cox_EM <- survival::coxph(Surv(time_event_EM, flg_event)~treat+X1+X2+X4+X5+X6+X8, data = df_EM)
# 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_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)
LM_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_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)
# 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_cox_MB <- survival::coxph(Surv(tstart,tstop,flg_event_MB)~flg_treat_MB+X1+X2+X4+X5+X6+X8, data = df3_MB2)
# 真の値の推定値の保存(Log HR)
LHR_true[i] <- gamma_0
LHR_naive[i] <- res_cox_naive$coefficients[1]
LHR_EM[i] <- res_cox_EM$coefficients[1]
LHR_LM_1yr[i] <- res_cox_LM_1yr$coefficients[1]
LHR_LM_2yr[i] <- res_cox_LM_2yr$coefficients[1]
LHR_MB[i] <- res_cox_MB$coefficients[1]
}
# 結果の整理#######################################################
df_LHR <- data.frame(
LHR_true = LHR_true,
LHR_naive = LHR_naive,
LHR_EM = LHR_EM,
LHR_LM_1yr = LHR_LM_1yr,
LHR_LM_2yr = LHR_LM_2yr,
LHR_MB = LHR_MB
)
df_LHR2 <-
df_LHR %>% mutate(
bias_naive = (LHR_naive-LHR_true),
bias_EM = (LHR_EM-LHR_true),
bias_LM_1yr = (LHR_LM_1yr-LHR_true),
bias_LM_2yr = (LHR_LM_2yr-LHR_true),
bias_MB = (LHR_MB-LHR_true),
SE_naive = (LHR_naive-LHR_true)^2,
SE_EM = (LHR_EM-LHR_true)^2,
SE_LM_1yr = (LHR_LM_1yr-LHR_true)^2,
SE_LM_2yr = (LHR_LM_2yr-LHR_true)^2,
SE_MB = (LHR_MB-LHR_true)^2
)
# ggplotで扱うため, long形式に
df_LHR2_long <-
df_LHR2 %>% rowid_to_column() %>%
pivot_longer(
cols = starts_with("LHR"),
names_to = "Method",
values_to = "LHR") %>%
select(rowid, Method, LHR) %>%
mutate(Method=str_remove(Method,"LHR_")) %>%
filter(Method != "true") %>%
mutate(
Method = factor(Method, levels = c("naive","EM","LM_1yr","LM_2yr","MB"))
)
# 推定結果をヒストグラム
p_est <- ggplot(data=df_LHR2_long, aes(x=LHR)) +
geom_histogram(colour="white", bins = 30) +
facet_grid(Method~.) +
geom_vline(aes(xintercept=gamma_0), colour="red") +
theme_bw() +
ggtitle(label = str_c("Histograms of Estimates : Log HR = ", round(gamma_0,2)))
# average of Bias
ave_bias <-
df_LHR2 %>% summarise(
across(starts_with("bias"), mean)
)
# MSE
MSE <-
df_LHR2 %>% summarise(
across(starts_with("SE"), mean)
)
list(
hist = p_est,
ave_bias = ave_bias,
MSE = MSE)
}シミュレーション回数を1000回, サンプルサイズを1000とし(論文では10000だが速度の都合で1000としている), ハザード比を0.5, 0.7, 1.0, 1.3, 1.5 と変更して実験.
nsim <- 1000 # number of simulations
N <- 1000 # sample size(論文では10000だが速度の問題で1000としている)
# Hazard Ratio(HR)を 0.5, 0.7, 1.0, 1.3, 1.5 変更
set.seed(2022)res_HR_0.5 <- sim_Mi2016(nsim=nsim, N=N, HR_0=0.5)
res_HR_0.7 <- sim_Mi2016(nsim=nsim, N=N, HR_0=0.7)
res_HR_1.0 <- sim_Mi2016(nsim=nsim, N=N, HR_0=1.0)
res_HR_1.3 <- sim_Mi2016(nsim=nsim, N=N, HR_0=1.3)
res_HR_1.5 <- sim_Mi2016(nsim=nsim, N=N, HR_0=1.5)真のハザード比が0.5のときの各手法による推定値のヒストグラム.
res_HR_0.5$hist
真のハザード比が0.7のときの各手法による推定値のヒストグラム.
res_HR_0.7$hist
真のハザード比が1.0のときの各手法による推定値のヒストグラム.
res_HR_1.0$hist
真のハザード比が1.3のときの各手法による推定値のヒストグラム.
res_HR_1.3$hist
真のハザード比が1.5のときの各手法による推定値のヒストグラム.
res_HR_1.5$hist
df_ave_bias <-
rbind.data.frame(
res_HR_0.5$ave_bias,
res_HR_0.7$ave_bias,
res_HR_1.0$ave_bias,
res_HR_1.3$ave_bias,
res_HR_1.5$ave_bias
) %>% mutate(
HR = c(0.5, 0.7, 1.0, 1.3, 1.5)
)
df_ave_bias_long <-
df_ave_bias %>%
pivot_longer(
cols = starts_with("bias"),
names_to = "Method",
values_to = "bias"
) %>% mutate(Method=str_remove(Method,"bias_")) %>%
mutate(
Method = factor(Method, levels = c("naive","EM","LM_1yr","LM_2yr","MB"))
)真のハザード比を変えていったときの,対数ハザード比の真の値と各手法の推定値とのズレの平均値の折れ線グラフ.
ggplot(data=df_ave_bias_long, aes(x=HR, y=bias, colour=Method)) +
geom_line() +
geom_point(aes(shape=Method)) +
ggtitle(label = "average bias of estimated log HR ")
df_ave_MSE <-
rbind.data.frame(
res_HR_0.5$MSE,
res_HR_0.7$MSE,
res_HR_1.0$MSE,
res_HR_1.3$MSE,
res_HR_1.5$MSE
) %>% mutate(
HR = c(0.5, 0.7, 1.0, 1.3, 1.5)
)
df_ave_MSE_long <-
df_ave_MSE %>%
pivot_longer(
cols = starts_with("SE"),
names_to = "Method",
values_to = "MSE"
) %>% mutate(Method=str_remove(Method,"SE_")) %>%
mutate(
Method = factor(Method, levels = c("naive","EM","LM_1yr","LM_2yr","MB"))
)真のハザード比を変えていったときの,対数ハザード比の真の値と各手法の推定値とのズレの二乗値の平均値(MSE)の折れ線グラフ.
ggplot(data=df_ave_MSE_long, aes(x=HR, y=MSE, colour=Method)) +
geom_line() +
geom_point(aes(shape=Method)) +
ggtitle(label = "MSE of estimated log HR ")
Mi et al. 2016のβがtreatmentもeventも1のもののみ, 今回扱っている. コードの関数の引数にβを入れて変更できるようにすれば他のβの場合も扱えるが, 面倒なのとそんなに個人的に重要ではなさそうなのでここではやらない. ランドマーク解析法とMantel-Byarの方法については, 結果は似たようになっているが, Exclusion Methodについては, 少しズレているような気がする. 今のところ, どこが間違っているのか(それとも間違っていないのか)分かっていない. HRが大きくなればなるほど, biasは減っている. HRが大きいということはtreatmentを受けたらすぐにeventが起きるということになる. Exclusion Methodがどうバイアスを引き起こしているか分からないが, treatmentの時間が短い分だけbiasが生じないようになっているのかもしれない. naiveな手法はHRが1から離れるほどbiasが大きくなっている. 正しくシミュレーションできているのか完全な自信はないので, 今後もう少し学びながら, 検討していくこととする. Naive MethodとExclusion Methodについて修正した結果, 論文のシミュレーションの結果図を全体的に良い感じに再現できたと思われる. 良かった.
ランドマーク解析も頑張っているとはいえ, やはりMantel-Byarの方法がHRに関わらず良い性能を示している. Exclusion MethodはNaive Methodよりはマシだが, バイアスが生じていることが分かる. Naive Methodは論外的にズレている. おそらく正しくシミュレーションをreplicateできていると思う.
【学習予定】
今回は以下の4つの文献を参考にした.
-
Mi et al. 2016
“Use of the landmark method to address immortal person-time bias in comparative effectiveness research: a simulation study”
https://pubmed.ncbi.nlm.nih.gov/27350312/
治療の定義の都合でアウトカムが起きえないフォローアップの期間によって生じるimmortal-time biasがあるが, バイアスを回避する解析方法として, Mantel-Byarの方法とlandmarkの方法がある.
上の2つの方法とnaiveな方法とexclusionの方法の4つをシミュレーションで評価する.
論文ではSASで書かれていてコードもない. これをRで再現実験する. -
Leemis et al. 1990
“Variate generation for accelerated life and proportional hazards models with time dependent covariates” https://www.sciencedirect.com/science/article/abs/pii/0167715290900529
時間依存性共変量付きの比例ハザードモデルと加速故障モデルの乱数生成方法を解説.
一様乱数uを-log(1-u)に変換して, これが0~, ハザード関数の逆関数で時間tに対応させて乱数を生成. -
Weberpals et al. 2018
“Comparative performance of a modified landmark approach when no time of treatment data are available within oncological databases: exemplary cohort study among resected pancreatic cancer patients” https://pubmed.ncbi.nlm.nih.gov/30214315/
ランドマーク解析やMantel-Byar methodのSASのコードが掲載されている. -
Gleiss 2018
“An unjustified benefit: immortal time bias in the analysis of time-dependent events”
https://onlinelibrary.wiley.com/doi/10.1111/tri.13081
immortal time biasについて学ぶのに参考になった.
学んでいる中で, 下の最近のレビューを発見した. Mi et al. 2016では比較されていない手法も上がっている.
-
Wang et al. 2022
“Statistical Methods for Accommodating Immortal Time: A Selective Review and Comparison” https://arxiv.org/abs/2202.02369 2022年のpre printのレビュー. target trial emulationへのclone-censor-weightの3 stepsのアプローチも掲載されている.
この手法は下の文献で平易に解説されている. -
Hernán 2018
“How to estimate the effect of treatment duration on survival outcomes using observational data”
https://pubmed.ncbi.nlm.nih.gov/29419381/
一通り目を通してみて, この手法のやり方はなんとなく分かったけど, なぜこれが上手くいくのか全く分かっていない.
これを今回のようにシミュレーションデータ作って, 試してみたい.
【追記：2022/07/11】
- Putter and van Houwelingen2017
“Understanding Landmarking and Its Relation with Time-Dependent Cox Regression.”
ランドマーク解析とCox回帰の関係性について書いている論文. 理解を促しそうなので, 追加しておく.
