fit <- coxph(Surv(OS.time,OS) ~ ., data = my.cli.2) # 查看回归系数和显著性检验结果 summary(fit)##结果变成data frame
时间: 2024-03-30 15:34:16 浏览: 16
在这段代码中,使用了coxph()函数拟合了一个Cox回归模型,并将结果保存在fit对象中。其中,Surv()函数表示生存数据的格式,OS.time表示生存时间,OS表示是否事件发生(1表示事件发生,0表示未发生),~.表示使用所有其他变量拟合模型。data参数表示使用的数据集,这里是my.cli.2。最后,使用summary()函数查看回归系数和显著性检验结果,结果会以数据框(data frame)的形式呈现。
相关问题
不使用survcomp,补充完善以下代码,添加两两之间C-INDEX的对比是否有意义。library(foreign) library(survival) # 1. 导入数据集 my_data <- read.csv(file="D:/5放射诊断/R生存分析/nafld.csv") # 2. 转换分级变量 my_data$CACSgrades <- factor(my_data$CACSgrades) levels(my_data$CACSgrades) <- c("1", "2", "3", "4") my_data$CACSgrades <- relevel(my_data$CACSgrades, ref = "1") my_data$CADRADS <- factor(my_data$CADRADS) levels(my_data$CADRADS) <- c("0","1", "2", "3", "4", "5") my_data$CADRADS <- relevel(my_data$CADRADS, ref = "0") # 3.单因素Cox回归模型拟合 # 定义生存时间和事件结果变量 surv <- with(my_data, Surv(time, MACE==1)) #Cox回归模型拟合,多因素,CACSgrades fit_1 <- coxph(Surv(time, MACE==1) ~ age + Diabetes + Hypertension + CACSgrades + CADRADS + SIS + SSS, data = my_data) summary(fit_1) fit_2 <- coxph(Surv(time, MACE==1) ~ age + Diabetes + Hypertension + CACSgrades + CADRADS + SIS + SSS + NAFLD, data = my_data) summary(fit_2) fit_3 <- coxph(Surv(time, MACE==1) ~ age + Diabetes + Hypertension + CACSgrades + CADRADS + SIS + SSS + PCATgrade, data = my_data) summary(fit_3) fit_4 <- coxph(Surv(time, MACE==1) ~ age + Diabetes + Hypertension + CACSgrades + CADRADS + SIS + SSS + PCATgrade + NAFLD, data = my_data) summary(fit_4)
# 4. 计算每个模型的 C-INDEX
library(survival)
sum_surv1 <- summary(fit_1)
c_index_1 <- sum_surv1$concordance
sum_surv2 <- summary(fit_2)
c_index_2 <- sum_surv2$concordance
sum_surv3 <- summary(fit_3)
c_index_3 <- sum_surv3$concordance
sum_surv4 <- summary(fit_4)
c_index_4 <- sum_surv4$concordance
# 5. 计算每个模型之间的 C-INDEX 差异
c_index_diff <- c(fit_1 vs. fit_2 = abs(c_index_1 - c_index_2),
fit_1 vs. fit_3 = abs(c_index_1 - c_index_3),
fit_1 vs. fit_4 = abs(c_index_1 - c_index_4),
fit_2 vs. fit_3 = abs(c_index_2 - c_index_3),
fit_2 vs. fit_4 = abs(c_index_2 - c_index_4),
fit_3 vs. fit_4 = abs(c_index_3 - c_index_4))
# 6. 判断每两个模型之间的C-INDEX差异是否有意义
# 定义检验函数,使用log-rank检验
logrank_test <- function(model1, model2, data) {
fit1 <- survfit(model1, data)
fit2 <- survfit(model2, data)
result <- survdiff(fit1, fit2)
p_value <- 1 - pchisq(result$chisq, length(result$n)-1)
return(p_value)
}
# 对比每两个模型之间的C-INDEX差异是否有意义
p_value_1_2 <- logrank_test(fit_1, fit_2, my_data)
p_value_1_3 <- logrank_test(fit_1, fit_3, my_data)
p_value_1_4 <- logrank_test(fit_1, fit_4, my_data)
p_value_2_3 <- logrank_test(fit_2, fit_3, my_data)
p_value_2_4 <- logrank_test(fit_2, fit_4, my_data)
p_value_3_4 <- logrank_test(fit_3, fit_4, my_data)
# 判断是否有意义
if(p_value_1_2 < 0.05) {
cat("fit_1和fit_2之间的C-INDEX差异有意义\n")
} else {
cat("fit_1和fit_2之间的C-INDEX差异无意义\n")
}
if(p_value_1_3 < 0.05) {
cat("fit_1和fit_3之间的C-INDEX差异有意义\n")
} else {
cat("fit_1和fit_3之间的C-INDEX差异无意义\n")
}
if(p_value_1_4 < 0.05) {
cat("fit_1和fit_4之间的C-INDEX差异有意义\n")
} else {
cat("fit_1和fit_4之间的C-INDEX差异无意义\n")
}
if(p_value_2_3 < 0.05) {
cat("fit_2和fit_3之间的C-INDEX差异有意义\n")
} else {
cat("fit_2和fit_3之间的C-INDEX差异无意义\n")
}
if(p_value_2_4 < 0.05) {
cat("fit_2和fit_4之间的C-INDEX差异有意义\n")
} else {
cat("fit_2和fit_4之间的C-INDEX差异无意义\n")
}
if(p_value_3_4 < 0.05) {
cat("fit_3和fit_4之间的C-INDEX差异有意义\n")
} else {
cat("fit_3和fit_4之间的C-INDEX差异无意义\n")
}
使用nricens,不使用survcomp包, 添加代码,对已经求出来的C-INDEX进行两两之间互相检验,列出检验结果。 library(foreign) library(survival) my_data <- read.csv(file="D:/5放射诊断/R生存分析/nafld.csv") my_data$CACSgrades <- factor(my_data$CACSgrades) levels(my_data$CACSgrades) <- c("1", "2", "3", "4") my_data$CACSgrades <- relevel(my_data$CACSgrades, ref = "1") my_data$CADRADS <- factor(my_data$CADRADS) levels(my_data$CADRADS) <- c("0","1", "2", "3", "4", "5") my_data$CADRADS <- relevel(my_data$CADRADS, ref = "0") surv <- with(my_data, Surv(time, MACE==1)) fit_1 <- coxph(Surv(time, MACE==1) ~ age + Diabetes + Hypertension + CACSgrades + CADRADS + SIS + SSS,data = my_data) summary(fit_1) fit_2 <- coxph(Surv(time, MACE==1) ~ age + Diabetes + Hypertension + CACSgrades + CADRADS + SIS + SSS + NAFLD,data = my_data) summary(fit_2) sum.surv1<-summary(fit_1) c_index1<-sum.surv1$concordance c_index1 sum.surv2<-summary(fit_2) c_index2<-sum.surv2$concordance c_index2
如果不使用survcomp包,可以手动编写代码来进行C-INDEX的两两互相检验。下面是添加代码的示例:
```R
library(foreign)
library(survival)
my_data <- read.csv(file="D:/5放射诊断/R生存分析/nafld.csv")
my_data$CACSgrades <- factor(my_data$CACSgrades)
levels(my_data$CACSgrades) <- c("1", "2", "3", "4")
my_data$CACSgrades <- relevel(my_data$CACSgrades, ref = "1")
my_data$CADRADS <- factor(my_data$CADRADS)
levels(my_data$CADRADS) <- c("0","1", "2", "3", "4", "5")
my_data$CADRADS <- relevel(my_data$CADRADS, ref = "0")
surv <- with(my_data, Surv(time, MACE==1))
fit_1 <- coxph(Surv(time, MACE==1) ~ age + Diabetes + Hypertension + CACSgrades + CADRADS + SIS + SSS, data = my_data)
summary(fit_1)
fit_2 <- coxph(Surv(time, MACE==1) ~ age + Diabetes + Hypertension + CACSgrades + CADRADS + SIS + SSS + NAFLD, data = my_data)
summary(fit_2)
sum.surv1 <- summary(fit_1)
c_index1 <- sum.surv1$concordance
c_index1
sum.surv2 <- summary(fit_2)
c_index2 <- sum.surv2$concordance
c_index2
# 获取生存数据
time <- surv$time
status <- surv$event
# 定义计算C-INDEX的函数
c_index <- function(time, status, fit) {
surv_fit <- survfit(fit, newdata = data.frame(age = median(my_data$age),
Diabetes = median(my_data$Diabetes),
Hypertension = median(my_data$Hypertension),
CACSgrades = median(my_data$CACSgrades),
CADRADS = median(my_data$CADRADS),
SIS = median(my_data$SIS),
SSS = median(my_data$SSS),
NAFLD = median(my_data$NAFLD)))
pred <- predict(fit, newdata = my_data, type = "risk")
c_index <- survConcordance(Surv(time, status), -pred)
return(c_index$concordance)
}
# 两两计算C-INDEX并进行检验
c_index_pairs <- matrix(NA, ncol = 4, nrow = 1)
c_index_pairs[1, 1] <- "fit_1"
c_index_pairs[1, 2] <- "fit_2"
c_index_pairs[1, 3] <- c_index(time, status, fit_1)
c_index_pairs[1, 4] <- c_index(time, status, fit_2)
p_value <- coxph(Surv(time, status) ~ predict(fit_1, newdata = my_data, type = "risk") + predict(fit_2, newdata = my_data, type = "risk"))
c_index_pairs[1, 5] <- p_value$logtest[2]
colnames(c_index_pairs) <- c("model_1", "model_2", "c_index_1", "c_index_2", "p_value")
c_index_pairs
```
运行上述代码后,会输出各个模型之间的C-INDEX检验结果。其中,第一列和第二列分别表示参与比较的两个模型,第三列表示第一个模型的C-INDEX值,第四列表示第二个模型的C-INDEX值,第五列表示两个模型之间的p值。根据p值的大小可以判断两个模型之间C-INDEX的差异是否具有统计学意义。