R语言survival包coxph函数拟合cox回归模型常用的科研论文参数的提取和合并
理论介绍：
Cox回归模型的基本形式如下：
H(t,x)=h0(t)exp(b’X)=h0(t)exp(b1x1+b2x2+...+bmxm)
其中：
X:指的就是自变量，一般有多个自变量
t：时间，也就是生存时间
H(t,X)：当生存时间为t时，自变量x时候的危险率。
我们通过会在队列研究中建立cox回归模型，然后根据估计的回归系数b1、b2、。。。bm取指数得到HR值，如x1变量的HR=eb1，x2变量的HR=eb2，即e的b1/b2次方。
我们在研究中，一般需要提取以下参数：
回归系数、回归系数标准误、检验统计量Z值、P值、HR值、HR值的置信区间
HR的和其置信区间合并，如1.34（1.22~1.58）
本文将基于已有数据：分别分析性别,年龄,城乡,文化程度对全因死亡、心血管死亡的影响，然后提取上述参数。
整个实践操作分为四步：
1.数据的整理：包括分类变量因子化
2.拟合模型，了解提取参数函数
3.建立函数，建立循环，提取两个模型参数并合并
4.结果导出
实践操作：
目的：分别分析性别,对全因死亡、心血管死亡的影响
自变量：性别
因变量：全因死亡、心血管死亡
协变量：年龄,城乡,文化程度
提取指标：回归系数、回归系数标准误、检验统计量Z值、P值、HR值、HR值的置信区间、
HR的和其置信区间的合并
因此，我们需要建立两个模型，分别是性别与全死因；性别与心血管死亡；然后提取两个模型中性别估计的参数
具体最终的结果如下：
代码实现
#########第一步：数据整理
.libPaths()#查看R包位置.libPaths("D:/Program File/R/R-4.3.2/library")#定义包安装位置setwd("E:/02学习/经验/03R语言图形绘制/17R语言survival包coxph函数拟合cox回归模型常用的科研论文参数的提取和合并")#设置工作空间getwd()#加载工作空间#读取要用的包library(survival)library(tableone) #提取HR(95%CI)和P
#第一步：导入数据,整理数据
#（1）导入数据data<-read.csv("data.csv",as.is = TRUE,header = T,sep = ",", fileEncoding='utf-8')#（2）将数据框的列名称转换为小写格式data <- setNames(data, tolower(colnames(data)))colnames(data)#（3）转因子a<-c("location","education")#填入需要转化的变量data[,a]<-lapply(data[,a],as.factor)#转因子
#第二步：建立模型，了解相关提取的函数
##cox回归
cox_fit1<- coxph(Surv(time, death_all) ~ gender+age+location+education,data=data)cox_fit2<- coxph(Surv(time, death_cvd) ~ gender+age+location+education,data=data)
##结果保存到a
a<-summary(cox_fit1)#下面（1, 3, 4)提取bata、标准误、z值
a$coefficients[, c(1, 3, 4)]#下面（1, 3, 4)表示提取HR值、HR的下限和HR的上限
a$conf.int[, c(1, 3, 4)]#这个是提取HR值和置信区间的合并，p值
ShowRegTable(cox_fit1,exp = TRUE,digits = 2,pDigits = 3,printToggle = TRUE, quote = FALSE, ciFun = confint)
#第三步
#（1）定义一个函数，用于提取主要暴露的参数估计值，这里的主要暴露是性别
# 定义一个函数，提取相关数据,整理到一个数据集中process_model <- function(cox_fit) {#结果保存为aa<-summary(cox_fit)##下面（1, 3, 4)提取bata、标准误、z值multi1 <- as.data.frame(round(a$coefficients[, c(1, 3, 4)], 2))#下面（1, 3, 4)表示提取HR值、HR的下限和HR的上限multi2 <- as.data.frame(round(a$conf.int[, c(1, 3, 4)], 2))#这个是提取HR值和置信区间的合并，p值multi3 <- ShowRegTable(cox_fit, exp = TRUE, digits = 2, pDigits = 3, printToggle = TRUE, quote = FALSE, ciFun = confint)##合并三个结果merged_result <- cbind(multi1, multi2,multi3)##将第一列命名为Characteristics；这是原来的变量merged_result <- tibble::rownames_to_column(merged_result, var = "Characteristics")##增加行，因为总共有9列，所以下面得加上NA需要有9列merged_result <- rbind(c("模型", NA, NA, NA, NA, NA,NA,NA,NA),#第1行merged_result[1,], #提取前1行，主要暴露是二分类，结果是1行c("对照", NA, NA, NA, NA, NA,NA,NA,NA))#增加最后1行return(merged_result)}
#（2）建立一个循环，将两个模型结果合并
final_result <- data.frame()for (model_name in c("cox_fit1", "cox_fit2")) {cox_fit <- get(model_name)current_result <- process_model(cox_fit)final_result <- rbind(final_result, current_result)}print(final_result)
#第四步：导出数据集
write.csv(final_result,file="结果.csv") 最重要的原始数据如下：供大家参考：
