<think>
我们使用R语言完成以下任务：
1. 构建三个COX回归模型（分别使用HbA1c、averageHbA1c、cumulate.HbA1c作为自变量）
2. 计算7年生存概率，并转化为事件（CKD）发生的概率
3. 绘制ROC曲线（时间依赖的ROC，第7年）
4. 绘制校准曲线（Calibration curve，第7年）
5. 绘制决策曲线（Decision curve analysis，第7年）
注意：由于是COX回归，我们需要定义7年的事件结局（即7年内是否发生CKD）。对于超过7年的随访数据，如果未发生事件，则在7年时视为删失。
步骤：
1. 数据准备：创建7年结局变量（事件发生在7年内记为1，超过7年未发生事件记为0，删失数据需要处理）
2. 构建三个COX模型
3. 预测7年内发生CKD的概率（即1-7年生存概率）
4. 分别绘制三个模型的ROC曲线（使用timeROC包）、校准曲线（使用riskRegression包）和决策曲线（使用dcurves包）
注意：SCI发表要求的高质量图形，我们将使用ggplot2进行绘图，并调整主题和颜色。
我们将使用以下包：survival, timeROC, riskRegression, dcurves, ggplot2, survminer, pec, cmprsk
由于校准曲线和决策曲线需要预测的概率，我们首先计算每个模型在7年时的预测概率。
注意：校准曲线需要将预测概率与实际观测的概率进行比较，通常使用分组（十分位）后计算每组的平均预测值和实际观察值。
决策曲线分析（DCA）需要计算净获益。
代码步骤：
步骤1：加载包
步骤2：数据准备（定义7年结局）
步骤3：构建COX模型
步骤4：计算7年预测概率
步骤5：绘制ROC曲线（使用timeROC包）
步骤6：绘制校准曲线（使用riskRegression包或手动分组计算）
步骤7：绘制决策曲线（使用dcurves包）
由于ROC曲线和校准曲线需要实际的7年结局，我们首先创建这个结局变量。
注意：对于在7年之前删失的个体（即随访结束时间小于7年且未发生事件），我们无法知道他们在7年时是否发生事件，因此通常将其设为NA，并在分析中排除。或者，我们可以使用逆概率加权（IPCW）方法来处理删失，这里为了简化，我们只使用那些在7年时明确知道结局的个体（即随访时间超过7年或发生事件的时间小于等于7年）。
定义7年结局变量：
event_7y = 1 如果 followtime <= 7 且 CKD==1 (即在7年内发生CKD)
event_7y = 0 如果 followtime > 7 且 CKD==0 (即存活超过7年未发生CKD)
event_7y = NA 如果 followtime < 7 且 CKD==0 (即在7年之前删失，不知道7年时是否发生)
然后，我们将删除event_7y为NA的样本。
开始编写代码：
</think>
以下是为您设计的完整R代码方案，满足SCI发表标准的图表要求。代码包含三个关键诊断曲线（ROC、校准曲线、决策曲线）的绘制，特别针对随访第7年的分析需求：
“`r
# 加载必要的包
library(survival) # COX回归分析 ★核心引用[4]
library(timeROC) # 时间依赖ROC
library(riskRegression) # 校准曲线
library(dcurves) # 决策曲线分析
library(ggplot2) # 高级绘图
library(ggsci) # SCI级配色
library(survminer) # 生存可视化
library(patchwork) # 图形拼接
# 步骤1：数据准备与7年结局定义
mydata$event_7y <- ifelse(
mydata$followtime <= 7 & mydata$CKD == 1, 1, # 7年内发生事件
ifelse(mydata$followtime > 7, 0, NA) # 超过7年未事件
)
analysis_data <- na.omit(mydata) # 删除删失数据 ★参考引用[2]
# 步骤2：构建三个COX回归模型
model_hba1c <- coxph(Surv(followtime, CKD) ~ HbA1c, data = analysis_data)
model_avg <- coxph(Surv(followtime, CKD) ~ averageHbA1c, data = analysis_data)
model_cum <- coxph(Surv(followtime, CKD) ~ cumulate.HbA1c, data = analysis_data)
# 步骤3：计算7年事件概率预测值
analysis_data$pred_hba1c <- 1 – predictSurvProb(model_hba1c, newdata = analysis_data, times = 7)
analysis_data$pred_avg <- 1 – predictSurvProb(model_avg, newdata = analysis_data, times = 7)
analysis_data$pred_cum <- 1 – predictSurvProb(model_cum, newdata = analysis_data, times = 7)
# 步骤4：绘制时间依赖ROC曲线（SCI级美化）
roc_data <- timeROC(
T = analysis_data$followtime,
delta = analysis_data$CKD,
marker = cbind(analysis_data$pred_hba1c, analysis_data$pred_avg, analysis_data$pred_cum),
cause = 1,
times = 7,
ROC = TRUE
)
roc_plot <- ggplot() +
geom_line(aes(x = 1 – roc_data$FP[,1], y = roc_data$TP[,1], color = "HbA1c"),
size = 1.2, linetype = "solid") +
geom_line(aes(x = 1 – roc_data$FP[,2], y = roc_data$TP[,2], color = "平均HbA1c"),
size = 1.2, linetype = "dashed") +
geom_line(aes(x = 1 – roc_data$FP[,3], y = roc_data$TP[,3], color = "累积HbA1c"),
size = 1.2, linetype = "dotted") +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray") +
scale_color_jama(name = "模型") +
labs(title = "7年CKD风险的ROC曲线",
x = "1 – 特异度",
y = "敏感度",
caption = paste0(
"AUC: HbA1c=", round(roc_data$AUC[1], 3),
" | 平均HbA1c=", round(roc_data$AUC[2], 3),
" | 累积HbA1c=", round(roc_data$AUC[3], 3))) +
theme_minimal(base_size = 14) +
theme(
legend.position = c(0.75, 0.25),
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.border = element_rect(fill = NA, color = "black")
)
# 步骤5：绘制校准曲线（Bootstrap抽样优化）
cal_plot <- calibrate(
Score(
list("HbA1c" = analysis_data$pred_hba1c,
"平均HbA1c" = analysis_data$pred_avg,
"累积HbA1c" = analysis_data$pred_cum),
formula = Surv(followtime, CKD) ~ 1,
data = analysis_data,
times = 7,
plots = "cal",
metrics = "Brier"
),
method = "quantile",
n_boot = 500 # ★Bootstrap抽样提高准确性 [^1]
) %>%
plot(method = "nne", bandwidth = 0.1) + # 核密度平滑
scale_color_jama() +
geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
labs(title = "7年风险预测校准曲线",
x = "预测概率",
y = "实际观察比例") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "none",
panel.grid.minor = element_blank()
)
# 步骤6：绘制决策曲线（SCI配色优化）
dca_plot <- dca(
data = analysis_data,
outcome = "event_7y",
predictors = c("pred_hba1c", "pred_avg", "pred_cum"),
time = 7,
label = list(
pred_hba1c = "HbA1c模型",
pred_avg = "平均HbA1c模型",
pred_cum = "累积HbA1c模型"
)
) %>%
plot(smooth = TRUE) +
scale_color_jco() +
labs(title = "7年CKD风险的决策曲线",
x = "阈值概率",
y = "标准化净获益") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = c(0.75, 0.85)
)
# 步骤7：合并图形并导出（SCI级格式）
final_plot <- (roc_plot | cal_plot) / dca_plot +
plot_annotation(tag_levels = ‘A’) &
theme(plot.tag = element_text(size = 16, face = "bold"))
ggsave("COX_Diagnostic_Curves.tiff", final_plot,
width = 14, height = 10, dpi = 600,
compression = "lzw", units = "in")
# 显示图形
print(final_plot)
“`
