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)險(xiǎn)表
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)險(xiǎn)表
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表
