后台粉丝希望我介绍一下能进行生存数据时点加权的ipwtm函数,今天简单分享一下。ipwtm函数属于ipw包,我们先导入包和数据。
library(ipw)
data(haartdat)
head(haartdat,6)
## patient tstart fuptime haartind event sex age cd4.sqrt endtime dropout
## 1 1 -100 0 0 0 1 22 23.83275 2900 0
## 2 1 0 100 0 0 1 22 25.59297 2900 0
## 3 1 100 200 0 0 1 22 23.47339 2900 0
## 4 1 200 300 0 0 1 22 24.16609 2900 0
## 5 1 300 400 0 0 1 22 23.23790 2900 0
## 6 1 400 500 0 0 1 22 24.85961 2900 0
这是一个HIV患者的生存数据,比较的是一种HAART治疗治疗方法对患者生存的影响。patient患者的ID,tstart开始随访的时间,fuptime随访终止的时间,haartind是否采用HAART治疗,event结局指标,是否死亡。sex性别,age随访时候的年龄,cd4.sqrt:CD4细胞的平方根。 ipwtm函数对生存数据加权很简单,其实就是一句话代码,exposure就是放入你的研究变量,numerator这里填入权重分子,如果按照Robins的算法这里填入1,属于不稳定权重,denominator这里属于权重分母,填入我们的协变量,它对权重的计算有点类似ipwpoint函数,不过更复杂,但是大概的原理是分子这个地方的变量生成一个模型,分母这个地方的变量生成一个模型,然后对预测值相除,当然过程中还有复杂的处理。tstart开始时间,timevar结束时间,data就是你的数据。type的类型有”first”, “cens” and “all”. “first”是权重由最低值进行转换。
temp <- ipwtm(
exposure = haartind,
family = "survival",
numerator = ~ sex + age,
denominator = ~ sex + age + cd4.sqrt,
id = patient,
tstart = tstart,
timevar = fuptime,
type = "first",
data = haartdat)
跑出结果后权重就是在temp中
我们可以把权重提取出来绘图
ipwplot(weights = temp$ipw.weights, timevar = haartdat$fuptime,
binwidth = 100, ylim = c(-1.5, 1.5), main = "Stabilized inverse probability weights")
这个图画法其实很简单,X就是我们的随访时间,Y轴就是权重取了个对数,设定一段宽度,每隔一段画一个箱线图。如果把ype设置为”cens”,进行的是逆概率删减权重(IPCW),我也不是很了解。
temp2 <- ipwtm(
exposure = dropout,
family = "survival",
numerator = ~ sex + age,
denominator = ~ sex + age + cd4.sqrt,
id = patient,
tstart = tstart,
timevar = fuptime,
type = "cens",
data = haartdat)
ipwplot(weights = temp2$ipw.weights, timevar = haartdat$fuptime,
binwidth = 100, ylim = c(-1.5, 1.5), main = "Stabilized inverse probability of censoring weights")
下面是加权和不加权对结果的影响
加权的
require(survival)
## 载入需要的程辑包:survival
summary(coxph(Surv(tstart, fuptime, event) ~ haartind + cluster(patient),
data = haartdat, weights = temp$ipw.weights*temp2$ipw.weights))
## Call:
## coxph(formula = Surv(tstart, fuptime, event) ~ haartind, data = haartdat,
## weights = temp$ipw.weights * temp2$ipw.weights, cluster = patient)
##
## n= 19175, number of events= 31
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## haartind -0.9363 0.3921 0.4299 0.4527 -2.068 0.0386 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## haartind 0.3921 2.551 0.1614 0.9521
##
## Concordance= 0.597 (se = 0.029 )
## Likelihood ratio test= 5.57 on 1 df, p=0.02
## Wald test = 4.28 on 1 df, p=0.04
## Score (logrank) test = 5.07 on 1 df, p=0.02, Robust = 4.83 p=0.03
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).
不加权的
summary(coxph(Surv(tstart, fuptime, event) ~ haartind, data = haartdat))
## Call:
## coxph(formula = Surv(tstart, fuptime, event) ~ haartind, data = haartdat)
##
## n= 19175, number of events= 31
##
## coef exp(coef) se(coef) z Pr(>|z|)
## haartind -0.5862 0.5565 0.4368 -1.342 0.18
##
## exp(coef) exp(-coef) lower .95 upper .95
## haartind 0.5565 1.797 0.2364 1.31
##
## Concordance= 0.569 (se = 0.032 )
## Likelihood ratio test= 1.97 on 1 df, p=0.2
## Wald test = 1.8 on 1 df, p=0.2
## Score (logrank) test = 1.85 on 1 df, p=0.2
加权的haartind是0.3921,不加权的haartind是0.5565,加权后影响变小。
绘制加权后的生存函数曲线
f1<-survfit(Surv(tstart, fuptime, event) ~ haartind + cluster(patient),
data = haartdat, weights = temp$ipw.weights*temp2$ipw.weights)
ggsurvplot(f1, data = haartdat)