Do optimal prognostic thresholds in continuous physiological variables really exist?
# https://scholar.google.com/scholar?hl=zh-TW&as_sdt=0%2C5&q=Do+optimal+prognostic+thresholds+in+continuous+physiological+variables+really+exist%3F+Analysis+of+origin+of+apparent+thresholds%2C+with+systematic+review+for+peak+oxygen+consumption%2C+ejection+fraction+and+BNP&btnG=#d=gs_qabs&t=1737163951359&u=%23p%3D171N1xvToIYJ
# 真實情況：風險因子每增加 1，死亡率線性增加 1%
# 結果：Kaplan-Meier、ROC 的閾值都類似於風險因子的中位數
# 設定文字編碼，確保輸出文字正常顯示
Sys.setlocale("LC_ALL", "C.UTF-8")
# 加載必要的套件
library(survival)
library(pROC)
# 設定模擬參數
set.seed(123)
n_patients <- 1500
risk_factor <- seq(0.01, 15.00, length.out = n_patients)
annual_mortality <- risk_factor / 100 # 風險因子（0.01, 15.00）每增加 1，死亡率線性增加 1%
# 計算風險因子中位數
median_risk_factor <- median(risk_factor)
# 模擬10年存活數據
simulate_survival <- function(annual_mortality, n_years = 10) {
survival_time <- numeric(length(annual_mortality))
event <- numeric(length(annual_mortality))
for (i in seq_along(annual_mortality)) {
for (year in 1:n_years) {
if (runif(1) < annual_mortality[i]) {
survival_time[i] <- year
event[i] <- 1
break
}
}
if (survival_time[i] == 0) {
survival_time[i] <- n_years
}
}
return(data.frame(survival_time, event))
}
# 生成模擬數據
sim_data <- simulate_survival(annual_mortality)
# Kaplan-Meier 分析函數
km_analysis <- function(data, risk_factor, thresholds) {
results <- data.frame(threshold = thresholds, chi_sq = numeric(length(thresholds)))
for (i in seq_along(thresholds)) {
data$group <- ifelse(risk_factor < thresholds[i], 0, 1)
if (length(unique(data$group)) > 1) {
chi_sq <- survdiff(Surv(survival_time, event) ~ group, data = data)$chisq
results$chi_sq[i] <- chi_sq
} else {
results$chi_sq[i] <- NA
}
}
return(results)
}
# 設定測試的風險因子閾值
thresholds <- seq(0.01, 15, by = 0.5)
# 執行Kaplan-Meier分析
km_results <- km_analysis(sim_data, risk_factor, thresholds)
valid_results <- km_results[!is.na(km_results$chi_sq), ]
best_threshold_km <- valid_results$threshold[which.max(valid_results$chi_sq)]
print(paste("最佳閾值 (Kaplan-Meier):", round(best_threshold_km, 2)))
# ROC 分析
roc_obj <- roc(sim_data$event, risk_factor)
best_threshold_roc <- coords(roc_obj, "best", ret = "threshold", transpose = FALSE)
print(paste("最佳閾值 (ROC):", round(best_threshold_roc, 2)))
# 繪製Kaplan-Meier分析結果
plot(valid_results$threshold, valid_results$chi_sq, type = "b",
xlab = "Threshold", ylab = "Chi-squared Value",
main = "Kaplan-Meier Analysis of Thresholds")
abline(v = best_threshold_km, col = "red", lty = 2, lwd = 2)
abline(v = median_risk_factor, col = "green", lty = 2, lwd = 2)
text(best_threshold_km, max(valid_results$chi_sq),
labels = paste0("KM Best: ", round(best_threshold_km, 2)),
pos = 4, col = "red")
text(median_risk_factor, max(valid_results$chi_sq) * 0.8,
labels = paste0("Median: ", round(median_risk_factor, 2)),
pos = 4, col = "green")
# 繪製ROC曲線
plot(roc_obj, main = "ROC Curve")
abline(v = best_threshold_roc, col = "blue", lty = 2, lwd = 2)
abline(v = median_risk_factor, col = "green", lty = 2, lwd = 2)
text(best_threshold_roc, 0.5,
labels = paste0("ROC Best: ", round(best_threshold_roc, 2)),
pos = 4, col = "blue")
text(median_risk_factor, 0.4,
labels = paste0("Median: ", round(median_risk_factor, 2)),
pos = 4, col = "green")
# 添加說明
legend("bottomright", legend = c("KM Best Threshold", "ROC Best Threshold", "Median Risk Factor"),
col = c("red", "blue", "green"), lty = 2, lwd = 2)
留言
張貼留言