R语言统计与绘图：COX回归模型怎么建？

install.packages("survival")
install.packages("survminer")
install.packages("forestplot")
install.packages("tableone")
install.packages("broom")
# 不需要的包可以不安装

library(survival)   # 生存分析需要
library(survminer)  # 作图使用ggplot画
library(forestplot)   # 画森林图
library(tableone)   
library(broom)   #  后面tidy函数需要
data("lung") 
head(lung)    # 预览 lung 数据集信息

fit <- survfit(Surv(time,status) ~ sex, data = lung);fit  # 以性别为例
survdiff(Surv(time,status) ~ sex,data = lung, rho = 0)  # 比较分析两组的生存时间分布是否不同
# rho = 0 表示使用long-rank检验或者Mantel-Haenszel 检验)

ggsurvplot(fit,
           data = lung,
           pval = TRUE,   # 显示P值,默认P值检验为 log-rank test.
           pval.coord = c(800,0.45),  #自定义P值在图中的坐标（x，y）,可自行调整
           conf.int = TRUE,    # 图中输出95%CI
           risk.table = TRUE,  # 输出risk table,“FALSE”为不显示；
           risk.table.height = 0.25,  # 调整 risk table的高度，默认为0.25
           risk.table.y.text = FALSE, # 不显示risk table分组文字，用颜色替代
           risk.table.y.text.col=TRUE,  # risk table文字颜色
           censor.shape = "+",  #可删除
           surv_median.line = "hv",   #输出中位生存期
           xlim = c(0,1100), break.x.by = 100,  # 自定义X轴范围和间距，自行调整
           ylim = c(0,1), break.y.by = 0.2, # 自定义y轴范围和间距，自行调整
           xlab = "Follow-up time(d)",  # 定义 X 轴标签
           font.x = c(14,"black"),  # x轴字体，颜色
           font.y = c(14,"black"),  # y轴字体，颜色
           font.tickslab = c(12, "plain", "black"),  # 标签字体，样式，颜色
           legend = c(0.85,0.75),  # 图例坐标位置（x，y），X,Y轴值都在0-1之间
           legend.title = "",  # 图例标题，在双引号中输入文字
           legend.labs = c("Male", "Female"))  #

res.cox1 <- coxph(Surv(time, status) ~ sex, data = lung)  # 以性别为例
summary(res.cox1)  # 生成完整的结果报告

covariates <- c("age", "sex",  "ph.karno", "ph.ecog", "wt.loss")  # 选中需要进行单因素分析的变量
univ_formulas <- sapply(covariates,
                        function(x) as.formula(paste('Surv(time, status)~', x)))
univ_formulas

univ_models <- lapply(univ_formulas, function(x){coxph(x, data = lung)})
univ_models    # 输出每个单因素COX回归分析结果

univ_results <- lapply(univ_models,
                       function(x){ 
                         x <- summary(x)
                         p.value<-signif(x$wald["pvalue"], digits=2)
                         wald.test<-signif(x$wald["test"], digits=2)
                         beta<-signif(x$coef[1], digits=2); #coeficient beta
                         HR <-signif(x$coef[2], digits=2);  #exp(beta)
                         HR.confint.lower <- signif(x$conf.int[,"lower .95"],2)
                         HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
                         HR <- paste0(HR, " (", 
                                      HR.confint.lower, "-", HR.confint.upper, ")")
                         res<-c(beta, HR, wald.test, p.value)
                         names(res)<-c("beta", "HR (95% CI for HR)",           "wald.test", "p.value")
                         return(res)
                         #return(exp(cbind(coef(x),confint(x))))
                       })
class(univ_results)
str(univ_results)

res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)

res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno, data =  lung)
summary(res.cox)   # 查看con回归模型结果

table1<-ShowRegTable(res.cox, exp = TRUE, digits = 2, pDigits = 3,
                      printToggle = TRUE, quote = FALSE, ciFun = confint)
table2<-tidy(res.cox)
table3<-cbind(table1,table2) #table1 table2是两个结果，现在合并成一个总表
write.csv(table3,file="table3.csv") # 数据导出

res.cox1 <- cox.zph(res.cox);res.cox1

temp<-cox.zph(res.cox)
ggcoxzph(temp)

cox3 <- coxph(Surv(time,status) ~ sex + age + ph.ecog + strata(pat.karno), data = lung); summary(cox3)  # 对pat.karno进行分层
res.cox2 <- cox.zph(cox3)   # 检验模型是不是等比例模型
res.cox2  # 查看结果

ggforest(res.cox,
         data = lung,
         main="hazard ratio",
         cpositions=c(0.02,0.22,0.4),
         fontsize=0.8,
         refLabel="reference",
         noDigits=2)