# 加载所需的包
library(survival)
library(timeROC)
library(readr)
library(dplyr)
# 设置随机种子以确保结果可重现
set.seed(123)
# 生成模拟数据
n <- 500 # 样本量
# 生成三个连续特征变量（假设这些特征与生存时间相关）
characteristic1 <- rnorm(n, mean = 0, sd = 1)
characteristic2 <- rnorm(n, mean = 2, sd = 1.5) # 调整均值和标准差
characteristic3 <- rnorm(n, mean = -1, sd = 0.8) # 调整均值和标准差
# 构建风险得分（线性组合）
risk_score <- 0.5 * characteristic1 + 1.2 * characteristic2 + 0.8 * characteristic3
# 生成生存时间（使用指数分布）
true_time <- rexp(n, rate = exp(-risk_score)) # 风险越高，生存时间越短
# 设置最大随访时间（删失时间）
censoring_time <- runif(n, min = 0, max = 36) # 随机删失时间，最大36个月
# 确定观察到的生存时间和状态
os_time <- pmin(true_time, censoring_time) # 观察到的时间是真实时间和删失时间的最小值
os_status <- as.numeric(true_time <= censoring_time) # 状态：1=事件发生，0=删失
# 创建数据框
data <- data.frame(
os_time = os_time,
os_status = os_status,
characteristic1 = characteristic1,
characteristic2 = characteristic2,
characteristic3 = characteristic3
)
# 查看数据基本情况
cat("数据基本统计信息：\n")
summary(data)
# 计算事件发生率和删失率
event_rate <- mean(data$os_status)
censoring_rate <- 1 - event_rate
cat(paste0("\n事件发生率: ", round(event_rate * 100, 2), "%\n"))
cat(paste0("删失率: ", round(censoring_rate * 100, 2), "%\n"))
# 构建Cox模型
f <- coxph(Surv(os_time, os_status) ~ characteristic1 + characteristic2 + characteristic3, data = data)
cat("\nCox模型结果:\n")
print(summary(f))
# 计算线性预测值
data$lp <- predict(f, newdata = data, type = "lp")
# 计算时间依赖ROC
cat("\n正在计算时间依赖ROC曲线...\n")
time_roc <- timeROC(
T = data$os_time, # 指定观察的生存时间
delta = data$os_status, # 生存结局
marker = data$lp, # 预测因子（线性预测值）
cause = 1, # 阳性结局指标值
weighting = "marginal", # 处理删失数据的方法
times = c(12, 18), # 计算12个月和18个月的ROC曲线
ROC = TRUE, # 保存敏感性和特异性
iid = TRUE # 计算置信区间
)
# 格式化AUC结果
twelve_months <- paste0("12个月 AUC (95%CI)=",
sprintf("%.3f", time_roc$AUC[1]), "(",
sprintf("%.3f", confint(time_roc, level = 0.95)$CI_AUC[1, 1]/100), "-",
sprintf("%.3f", confint(time_roc, level = 0.95)$CI_AUC[1, 2]/100), ")")
eighteen_months <- paste0("18个月 AUC (95%CI)=",
sprintf("%.3f", time_roc$AUC[2]), "(",
sprintf("%.3f", confint(time_roc, level = 0.95)$CI_AUC[2, 1]/100), "-",
sprintf("%.3f", confint(time_roc, level = 0.95)$CI_AUC[2, 2]/100), ")")
# 输出AUC结果
cat("\nAUC结果:\n")
cat(paste0(twelve_months, "\n"))
cat(paste0(eighteen_months, "\n"))
# 绘制ROC曲线
if (!is.null(time_roc)) {
cat("\n正在绘制ROC曲线...\n")
plot(title = "", time_roc, col = "DodgerBlue", time = 12, lty = 1, lwd = 2)
plot(time_roc, time = 18, lty = 1, lwd = 2, add = TRUE, col = "LightSeaGreen")
legend("bottomright", c(twelve_months, eighteen_months),
col = c("DodgerBlue", "LightSeaGreen"), lty = 1, lwd = 2)
} else {
warning("time_roc为空，无法绘制ROC曲线")
}
为什么这个代码，可以运行，只要换成我的代码， 就会出现下标越界和In max(abs(colMeans(temp2))) : max里所有的参数都不存在；返回-Inf，# 加载所需的包
library(survival)
library(timeROC)
library(readr)
library(dplyr)
# 读取数据
data <- read_csv("data/共病数据/合并_累积结果/中风/HDL合并结果/new_folder/calculated_merged_data_final2_累积.csv")
# 检查关键变量是否存在
required_vars <- c("time", "status", "TyG_2011", "TyG_WC_2011", "TyG_WHtR_2011")
missing_vars <- required_vars[!required_vars %in% names(data)]
if (length(missing_vars) > 0) {
stop(sprintf("数据中缺少变量: %s", paste(missing_vars, collapse = ", ")))
}
# 构建Cox模型
if (nrow(data) > 0) {
f <- coxph(Surv(time, status) ~ TyG_2011 + TyG_WC_2011 + TyG_WHtR_2011, data = data)
} else {
stop("数据为空，无法构建Cox模型")
}
# 计算线性预测值
data$lp <- predict(f, newdata = data, type = "lp")
# 计算时间依赖ROC
time_roc <- timeROC(
T = data$time, # 指定观察的生存时间
delta = data$status, # 生存结局
marker = data$lp, # 预测因子（线性预测值）
cause = 1, # 阳性结局指标值
weighting = "marginal", # 处理删失数据的方法
times = c(12, 18), # 计算12个月和18个月的ROC曲线
ROC = TRUE, # 保存敏感性和特异性
iid = TRUE # 计算置信区间
)
# 格式化AUC结果
ci_result <- confint(time_roc, level = 0.95)$CI_AUC
if (nrow(ci_result) >= 1 && ncol(ci_result) >= 1) {
twelve_months <- paste0("12个月 AUC (95%CI)=",
sprintf("%.3f", time_roc$AUC[1]), "(",
sprintf("%.3f", ci_result[1, 1]/100), "-",
sprintf("%.3f", ci_result[1, 2]/100), ")")
eighteen_months <- paste0("18个月 AUC (95%CI)=",
sprintf("%.3f", time_roc$AUC[2]), "(",
sprintf("%.3f", ci_result[2, 1]/100), "-",
sprintf("%.3f", ci_result[2, 2]/100), ")")
} else {
warning("置信区间计算结果结构异常，无法正常提取值")
twelve_months <- NA
eighteen_months <- NA
}
# 输出AUC结果
cat("\nAUC结果:\n")
cat(paste0(twelve_months, "\n"))
cat(paste0(eighteen_months, "\n"))
# 绘制ROC曲线
if (!is.null(time_roc)) {
cat("\n正在绘制ROC曲线...\n")
plot(title = "", time_roc, col = "DodgerBlue", time = 12, lty = 1, lwd = 2)
plot(time_roc, time = 18, lty = 1, lwd = 2, add = TRUE, col = "LightSeaGreen")
legend("bottomright", c(twelve_months, eighteen_months),
col = c("DodgerBlue", "LightSeaGreen"), lty = 1, lwd = 2)
} else {
warning("time_roc为空，无法绘制ROC曲线")
}
本文介绍FXC中有关有时间依赖属性的主档维护方法，包括如何正确生成和维护这些属性，并讨论了查询时如何指定时间来获取准确信息。
扫一扫
1023
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
