基于R语言逆概率加权(IPTW)并行生存曲线分析

逆概率加权法最早由 Ｈｏｒｖｉｔｚ和Ｔｈｏｍｐｓｏｎ提出，即对每个可观测的ｙi的概率取倒数，作为被观测的 yi 的权重，修正由缺失数据或有偏抽样带来的估计偏差.IPTW 是减少多组观察性数据间混杂偏倚的有效方法, 在处理多组间变量混杂偏倚中起到了重要作用。简单来说，就是把许多协变量和混杂因素打包成一个概率并进行加权，这样的话，我只用计算它的权重就可以了，方便了许多。


经我自己总结，做R语言逆概率加权(IPTW)需要以下几步：

1. 导入数据，整理好数据
2. 把需要比较的变量拿出来做结果变量，其他的变量做协变量，组成logistic回归方程，生成预测值
3. 通过预测值计算出拟概率权重
4. 建立COX回归方程，加入权重
   今天我们使用R语言来演示逆概率加权(IPTW)并行生存曲线分析，继续使用我们的乳腺癌数据（公众号回复：乳腺癌可以获得数据），可以做逆概率加权(IPTW)的R包很多，今天我们以survival包和RISCA包分别做演示
   首先导入包和导入数据

library(foreign)
library(RISCA)
library(survminer)
bc <- read.spss("E:/r/test/Breast cancer survival agec.sav",
                use.value.labels=F, to.data.frame=T)
bc <- na.omit(bc)

我们先来看看数据：
age表示年龄，pathsize表示病理肿瘤大小（厘米），lnpos表示腋窝淋巴结阳性，histgrad表示病理组织学等级，er表示雌激素受体状态，pr表示孕激素受体状态，status结局事件是否死亡，pathscat表示病理肿瘤大小类别（分组变量），ln_yesno表示是否有淋巴结肿大，time是生存时间，后面的agec是我们自己设定的，不用管它。

把分类变量转换成因子

bc$er<-as.factor(bc$er)
bc$pr<-as.factor(bc$pr)
bc$ln_yesno<-as.factor(bc$ln_yesno)
names(bc)

OK，整理好数据之后，假设我们想知道ln_yesno（表示是否有淋巴结肿大）对生存时间的影响。
为了好进行比较，我们先做下没有概率加权(IPTW)的COX回归模型,可以看到ln_yesno1为1.9851

fit=coxph(Surv(time,status) ~ ln_yesno+ age+er+pr+histgrad+pathsize,data=bc)
summary(fit)

绘制生存曲线图

fit1 <- survfit(Surv(time,status) ~ ln_yesno, data = bc)
plot(fit1, ylab="Confounder-adjusted survival",
     xlab="Time post-transplantation (years)", col=c(1,2), grid.lty=1)

也可以使用ggsurvplot函数绘图

ggsurvplot(fit1, data = bc)

进一步修饰一下，可以看出有淋巴转移的生存率明显要低

下面我们来进行逆概率加权
我们把ln_yesno作为结果变量， age、er、pr、histgrad、pathsize作为协变量组成logistic回归方程

pr<- glm(ln_yesno~ age+er+pr+histgrad+pathsize, data=bc,
           family=binomial(link = "logit"))

生成预测值

pr1<-pr$fitted.values

生成权重

W <- (bc$ln_yesno==1) * (1/pr1) + (bc$ln_yesno==0) * (1)/(1-pr1)

查看生成权重后模型, ln_yesno1为1.9525，和没有拟加权之前（1.9851 ）结果有一点改变，但是改变不大，总体说明我们的结果是稳定的，可靠的

fit.IPTW=coxph(Surv(time,status) ~ ln_yesno+ age+er+pr+histgrad+pathsize,
               data=bc,weights=W)
summary(fit.IPTW)

生成权重后就可以用RISCA 包的ipw.survival函数可以进行加权并绘图了，和没有加权之前有一点点改变，红线稍微上移了一点

fit.ipw<-ipw.survival(times=bc$time, failures=bc$status,
                       variable=bc$ln_yesno, weights=W)
plot(fit, ylab="Confounder-adjusted survival",
     xlab="Time post-transplantation (years)", col=c(1,2), grid.lty=1)

也可以使用ggsurvplot函数绘图

fit2.ipw <- survfit(Surv(time,status) ~ ln_yesno, data = bc,weights = W)
ggsurvplot(fit2.ipw, data = bc)

OK,看着麻烦，其实很简单，主要就那几个步骤。