主题与背景
本文主要介绍了如何在R语言中使用predict函数对已拟合的Cox比例风险模型进行新数据的预测和推理。Cox模型是一种常用的生存分析方法，用于评估多个因素对事件发生时间的影响。文章通过具体的代码示例展示了如何使用predict函数的不同参数来获取生存概率和风险比。
主要观点
导入survival包
首先需要导入R语言中的survival包，该包提供了实现Cox比例风险模型和其他生存分析方法的功能。
假设已有Cox模型
假设已经有一个Cox模型，并将其存储在cox_model对象中。这个模型是通过之前的数据拟合得到的。
创建新数据集
为了进行预测，需要创建一个新的数据集new_data。在这个例子中，新数据集包含两列：年龄（age）和性别（sex），分别对应三个新的个体。
使用predict函数进行预测
生存概率预测：
使用predict函数并设置type = "survival"来预测新数据的生存概率。这将返回每个新个体在特定时间点的生存概率。
风险比预测：
使用predict函数并设置type = "risk"来预测新数据的风险比。这将返回每个新个体的风险比，表示相对于参考水平的风险增加或减少的程度。
打印预测结果
最后，文章展示了如何打印出预测的生存概率和风险比，以便进一步分析和解释。
# 导入必要的包
library(survival)# 创建一个包含更多变量的数据集
set.seed(123) # 为了可重复性设置随机数种子
data <- data.frame(time = c(5, 8, 10, 4, 12, 7, 9, 6, 13, 11), # 生存时间status = c(1, 0, 1, 0, 1, 1, 0, 1, 0, 1), # 结局状态：1表示事件发生，0表示删失age = c(60, 65, 70, 55, 80, 72, 60, 68, 75, 65), # 年龄sex = c("Male", "Female", "Male", "Female", "Male", "Female", "Male", "Female", "Male", "Female"), # 性别smoking_status = c("Smoker", "Non-smoker", "Non-smoker", "Smoker", "Smoker", "Non-smoker", "Smoker", "Non-smoker", "Smoker", "Non-smoker"), # 吸烟状态treatment = c("DrugA", "DrugB", "DrugA", "DrugA", "DrugB", "DrugA", "DrugB", "DrugA", "DrugB", "DrugA") # 治疗类型
)# 将分类变量转化为因子变量
data$sex <- factor(data$sex, levels = c("Male", "Female"))
data$smoking_status <- factor(data$smoking_status, levels = c("Smoker", "Non-smoker"))
data$treatment <- factor(data$treatment, levels = c("DrugA", "DrugB"))# 拟合Cox回归模型
cox_model <- coxph(Surv(time, status) ~ age + sex + smoking_status + treatment, data = data)# 打印Cox模型结果
summary(cox_model)# 创建新的数据集进行预测
new_data <- data.frame(age = c(60, 65, 70),sex = c("Male", "Female", "Male"),smoking_status = c("Smoker", "Non-smoker", "Smoker"),treatment = c("DrugA", "DrugB", "DrugA")
)# 将新的数据集的分类变量转换为因子
new_data$sex <- factor(new_data$sex, levels = c("Male", "Female"))
new_data$smoking_status <- factor(new_data$smoking_status, levels = c("Smoker", "Non-smoker"))
new_data$treatment <- factor(new_data$treatment, levels = c("DrugA", "DrugB"))# 使用predict函数进行生存概率和风险预测
predicted_survival <- predict(cox_model, newdata = new_data, type = "survival")
predicted_hazard <- predict(cox_model, newdata = new_data, type = "risk")# 打印预测结果
cat("预测的生存概率：\n")
print(predicted_survival)cat("\n预测的风险比：\n")
print(predicted_hazard)
总结
本文详细介绍了在R语言中使用predict函数对Cox比例风险模型进行新数据预测的具体步骤。核心观点包括导入必要的包、准备新数据集、使用predict函数的不同参数（type = "survival" 和 type = "risk"）来进行生存概率和风险比的预测，以及如何输出和查看这些预测结果。通过这些步骤，用户可以有效地利用已有的Cox模型对新数据进行生存分析和风险评估。
cox 代码如下：
# 加载所需的库
library(rms)
library(timeROC)# 假设我们的随访时间变量是 time, 结局状态变量是 status，并且我们有两个预测变量 predictor1 和 predictor2如果是分类变量，需要设置为factor, 并加上文字标签，连续性变量不用管，例如：
train_data$predictor1<-factor(train_data$predictor1,levels = c(0,1),labels = c('No','Yes'))
train_data$predictor2<-factor(train_data$predictor2,levels = c(1,2,3),labels = c('Stage I','Stage II','Stage III'))# 设置模型公式
formula <- Surv(time,status) ~ predictor1 + predictor2# 首先，我们需要为我们的数据设置一个数据分布对象
# 这会帮助rms包更好地理解我们的数据结构
options(datadist=NULL)
ddist <- datadist(train_data)
options(datadist = 'ddist')# 现在我们可以在训练数据上拟合Cox回归模型
model <- rms::cph(formula, data = train_data, surv = T,x = TRUE, y = TRUE)# 在训练集上生成预测值
train_predvalue <- predict(model, newdata = train_data)# 在训练集上生成ROC对象roc_train <- timeROC::timeROC(T=train_data$time,delta=train_data$status,marker=train_predvalue,cause=1,weighting='marginal',times=c(365.25*1,365.25*3,365.25*5),iid=TRUE)# 求auc和可信区间
roc_train$AUC
confint(roc_train)$CI_AUC# 绘制曲线plot(roc_train,time=365.25*1,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')plot(roc_train,time=365.25*3,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')plot(roc_train,time=365.25*5,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')# 在内部验证集上生成预测值
validation_predvalue <- predict(model, newdata = validation_data)# 在训练集上生成ROC对象roc_validation <- timeROC::timeROC(T=validation_data$time,delta=validation_data$status,marker=validation_predvalue,cause=1,weighting='marginal',times=c(365.25*1,365.25*3,365.25*5),iid=TRUE)# 求auc和可信区间
roc_validation$AUC
confint(roc_validation)$CI_AUC# 绘制曲线plot(roc_validation,time=365.25*1,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')plot(roc_validation,time=365.25*3,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')plot(roc_validation,time=365.25*5,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')# 在外部验证集上生成预测值
external_predvalue <- predict(model, newdata = external_data)# 在训练集上生成ROC对象roc_external <- timeROC::timeROC(T=external_data$time,delta=external_data$status,marker=external_predvalue,cause=1,weighting='marginal',times=c(365.25*1,365.25*3,365.25*5),iid=TRUE)# 求auc和可信区间
roc_external$AUC
confint(roc_external)$CI_AUC# 绘制曲线plot(roc_external,time=365.25*1,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')plot(roc_external,time=365.25*3,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')plot(roc_external,time=365.25*5,lty=1,lwd=2,xlab='1-Specificities',ylab='Sensitivities',col = 'red',title='')
