这是我在 stackoverflow 上的第一篇文章。我试图根据指南提出问题。
我目前正在研究一个边际结构模型。我已经阅读了 Hernan 关于这个主题的大部分文章,并阅读了他的书What if和提供的 R 代码。但是,我发现的所有示例都相当简单。例如,患者在随访开始时接受或不接受治疗 A 并在整个期间保持该治疗。我发现很难将我学到的知识应用到我的数据集上。
我的数据集是关于事件发生时间的数据。在患者的随访期间,他们可以接受治疗A。一些患者在随访开始时已经接受治疗A,而其他患者在随访期间较晚开始或根本不接受治疗。此外,正在接受治疗A的患者可以在他们的随访期间停止治疗并再次重新开始治疗,即随时间变化的治疗。
在这里,我提供了一个示例数据集,其中包含事件时间数据的粗略 MSM 代码。我使用合并逻辑回归代替 Cox PH 回归,其优势比类似于 Cox 模型的风险比。为简单起见,我没有计算审查权重。
set.seed(1948)
idvar<-1:100 #100 subjects
v1<-sample(c(0,1), replace=TRUE, size=length(idvar)) #One baseline variable
time<-sample(1:100, replace=TRUE,size=length(idvar)) #Follow-up time
event<-sample(c(0,1), replace=TRUE, size=length(idvar)) #End of follow-up event. If 0 patients either reached administrative censoring or a competing event
df<-data.frame(idvar,v1,time,event) #Create dataframe
library(splitstackshape)
df.expand<-expandRows(df,"time", drop = F) #Create dataframe with 1 row per 1 unit of time per subject
df.expand$time_1 <- sequence(rle(df.expand$idvar)$lengths)-1 #Calculates times in expanded dataframe
df.expand$event_true <- as.factor(ifelse(df.expand$time_1==df.expand$time-1 & #Only event at last follow-up
df.expand$event==1, 1, 0))
df.expand$treatment<-sample(c(0,1), replace=TRUE,size=nrow(df.expand)) #Time-varying treatment
df.expand$v2<-sample(c(0,1), replace=TRUE,size=nrow(df.expand)) #Time-varying covariate
#Calculate weights. Need help here
pred_dnom = predict(glm(treatment~factor(v1) + factor(v2), #Time-varying covariates only in denominator
data = df.expand,
family = binomial(link = "logit")), type = "response")
pred_num = predict(glm(treatment~factor(v1),
data = df.expand,
family = binomial(link = "logit")), type = "response")
df.expand$sw <- ifelse(df.expand$treatment == 0, ((1 - pred_num) / (1 - pred_dnom)),
(pred_num / pred_dnom))
##Normally I would also calculate censoring weights and multiply those with the treatment weights##
fit1<-glm(event_true==1~treatment+factor(time_1),
data = df.expand, weights = sw,
family = binomial(link = "logit"))
然而,在上面的例子中,权重不考虑以前的治疗史和随时间变化的混杂因素治疗史。所以我的问题是,如何将这些历史纳入时变治疗。Fewell等人的这篇文章。描述了针对时变治疗但患者无法退出并重新开始治疗的方法。我引用:
为了估计每个受试者到每个月的完整治疗史的概率(3中的分母),我们将每个月观察到的治疗的估计概率乘以时间累积。每个受试者的第一个估计概率保持不变。对于所有其他情况,当前时间点的估计概率乘以前一个时间点的估计概率。[Fewell Z,Hernán MA,Wolfe F,Tilling K,Choi H,Sterne JAC。使用边际结构模型控制时间相关的混杂。统计杂志。2004;4(4):402-420。doi:10.1177/1536867X0400400403]
从他们文章的代码中可以看出,一旦患者开始治疗,他们的权重就设置为 1。我知道很多人用来计算权重的另一种方法是ipw 包。但是,当患者开始治疗时,此软件包还将权重设置为 1,之后权重也保持为 1。
ipw 包的一位作者也为Grafféo 等人的文章做出了贡献。它扩展了 ipw 包以允许随时间变化的权重。他们通过计算在不接受治疗时开始治疗的权重和在接受治疗时停止治疗的权重来做到这一点。(代码可作为文章中的支持信息)但是,他们没有考虑治疗的历史/随时间变化的混杂因素。
更准确地说,我们没有考虑混杂治疗暴露的整个历史,而只是认为在时间 t 的感兴趣的暴露仅取决于在时间 t 给出的混杂治疗。然后,不是计算直到第一次修改感兴趣的曝光时的权重,而是计算所有时间点的权重,即,计算感兴趣状态处理的所有变化和所有事件时间的权重。[Grafféo, N, Latouche, A, Geskus, RB, Chevret, S. 使用治疗权重的反概率对时变暴露进行建模。生物识别杂志。2018; 60: 323–332. https://doi.org/10.1002/bimj.201600223]
他们为此提出以下建议
我们没有考虑混杂因素的历史。对于感兴趣的暴露和混杂因素之间的关系中的每个转换,可以使用不同的模型来解决这个问题。需要进一步的工作来研究这种建模选项。最简单的方法是将先前暴露的计数器作为协变量包含在治疗模型中。[Grafféo, N, Latouche, A, Geskus, RB, Chevret, S. 使用治疗权重的反概率对时变暴露进行建模。生物识别杂志。2018; 60: 323–332. https://doi.org/10.1002/bimj.201600223]
因此,Fewell 的方法确实结合了治疗和随时间变化的变量历史,但不允许在开始治疗后改变体重。Grafféo 的方法确实允许权重随时间变化,但不包含处理和随时间变化的变量历史。所以我真正想要的是结合这些方法,但我绝对不知道这怎么可能。
我希望我的问题很清楚(并且问得正确)并且有人有建议。谢谢!