library(survminer) # 加載包,沒有的包就按安裝
library(survival) 
library(tableone)
library(survey)
library(MatchIt)
library(reportReg)
library(foreign)

testdata<-read.csv2("seerdev.csv",header = T,sep = ",")#讀取外部csv格式數(shù)據(jù)并將數(shù)據(jù)賦值給testdata

testdata$age<-as.numeric(testdata$age)#把a(bǔ)ge變?yōu)閿?shù)值型變量

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

str(testdata)#查看數(shù)據(jù)集結(jié)構(gòu)

#下面首先展示未調(diào)整的生存曲線

fit <- survfit(Surv(OS,status) ~ sex,  # 創(chuàng)建生存對象 
               data = testdata) # 數(shù)據(jù)集來源
fit # 查看擬合曲線信息

ggsurvplot(fit, # 創(chuàng)建的擬合對象
           data = testdata,  # 指定變量數(shù)據(jù)來源
           conf.int = FALSE, # 顯示置信區(qū)間
           pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存時間線
           risk.table = TRUE, # 添加風(fēng)險表
           xlab = "Follow up time(d)", # 指定x軸標(biāo)簽
           legend = c(0.8,0.75), # 指定圖例位置
           legend.title = "", # 設(shè)置圖例標(biāo)題
           legend.labs = c("Male", "Female"), # 指定圖例分組標(biāo)簽
           break.x.by = 10)  # 設(shè)置x軸刻度間距





#后續(xù)進(jìn)行PSM及IPTW時以sex為分組變量，age和BMIteam為待調(diào)整變量計(jì)算PSM及進(jìn)行IPTW

attach(testdata)#加載數(shù)據(jù)集到環(huán)境中

vars<-c("age","BMIteam")#待調(diào)整變量組成向量并定義為vars

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

testdata$ps=predict(psModel,type="response")#計(jì)算傾向性評分并在數(shù)據(jù)集內(nèi)添加一列為PS的列，內(nèi)容為評分

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,# 創(chuàng)建生存對象 
                         data = testdata) # 數(shù)據(jù)集來源
summary(fit.IPTW)

ggsurvplot(fit.IPTW, # 創(chuàng)建的擬合對象
           data = testdata,  # 指定變量數(shù)據(jù)來源
           conf.int = FALSE, # 顯示置信區(qū)間
           pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存時間線
           risk.table = TRUE, # 添加風(fēng)險表
           xlab = "Follow up time(d)", # 指定x軸標(biāo)簽
           legend = c(0.8,0.75), # 指定圖例位置
           legend.title = "", # 設(shè)置圖例標(biāo)題
           legend.labs = c("Male", "Female"), # 指定圖例分組標(biāo)簽
           break.x.by = 10)  # 設(shè)置x軸刻度間距

#以下分別進(jìn)行單因素及多因素cox回歸分析

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

單因素HR表