R语言超详细逆概率加权(IPTW)后生存曲线做法

library(survminer) # 加载包,没有的包就按安装
library(survival) 
library(tableone)
library(survey)
library(MatchIt)
library(reportReg)
library(foreign)

testdata<-read.csv2("seerdev.csv",header = T,sep = ",")#读取外部csv格式数据并将数据赋值给testdata

testdata$age<-as.numeric(testdata$age)#把age变为数值型变量

testdata$sex<-factor(testdata$sex,labels=c("male","female"))#把sex变量因子化

str(testdata)#查看数据集结构

#下面首先展示未调整的生存曲线

fit <- survfit(Surv(OS,status) ~ sex,  # 创建生存对象 
               data = testdata) # 数据集来源
fit # 查看拟合曲线信息

ggsurvplot(fit, # 创建的拟合对象
           data = testdata,  # 指定变量数据来源
           conf.int = FALSE, # 显示置信区间
           pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存时间线
           risk.table = TRUE, # 添加风险表
           xlab = "Follow up time(d)", # 指定x轴标签
           legend = c(0.8,0.75), # 指定图例位置
           legend.title = "", # 设置图例标题
           legend.labs = c("Male", "Female"), # 指定图例分组标签
           break.x.by = 10)  # 设置x轴刻度间距





#后续进行PSM及IPTW时以sex为分组变量，age和BMIteam为待调整变量计算PSM及进行IPTW

attach(testdata)#加载数据集到环境中

vars<-c("age","BMIteam")#待调整变量组成向量并定义为vars

psModel<-glm(team~age+BMIteam,family=binomial(link="logit"),data=testdata) 

testdata$ps=predict(psModel,type="response")#计算倾向性评分并在数据集内添加一列为PS的列，内容为评分

head(testdata$ps)#展示前6个患者评分

testdata$IPTW<-ifelse(testdata$sex=="female",1/testdata$ps,1/(1-testdata$ps))

fit.IPTW<- survfit(Surv(OS,status) ~ sex, 
                         weights=testdata$IPTW,# 创建生存对象 
                         data = testdata) # 数据集来源
summary(fit.IPTW)

ggsurvplot(fit.IPTW, # 创建的拟合对象
           data = testdata,  # 指定变量数据来源
           conf.int = FALSE, # 显示置信区间
           pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存时间线
           risk.table = TRUE, # 添加风险表
           xlab = "Follow up time(d)", # 指定x轴标签
           legend = c(0.8,0.75), # 指定图例位置
           legend.title = "", # 设置图例标题
           legend.labs = c("Male", "Female"), # 指定图例分组标签
           break.x.by = 10)  # 设置x轴刻度间距

#以下分别进行单因素及多因素cox回归分析

model=coxph(Surv(OS,status)~sex,data=testdata)
model.IPTW=coxph(Surv(OS,status)~sex,data=testdata,weights=testdata$IPTW)