# -*- coding: utf-8 -*-
# 课题二：投稿前敏感性分析（回应 nature-reviewer 审稿意见）
suppressPackageStartupMessages({library(dplyr); library(readr)})

DATA <- "D:/龙虾工作空间/2026-08-07-17-00-35/监测数据/选题二_CKM危险因素聚集/课题二_分析数据.csv"
dat <- read_csv(DATA, show_col_types = FALSE)
dat$sex <- ifelse(dat$性别 == "男", "Male", "Female")
dat$sex <- factor(dat$sex, levels = c("Male", "Female"))
dat$年龄 <- as.numeric(dat$年龄)

# ---------- 结局（对齐选题一） ----------
dat$qtc_prol <- ifelse(dat$sex == "Male" & dat$ECG_QTc_ms >= 450, 1,
                ifelse(dat$sex == "Female" & dat$ECG_QTc_ms >= 460, 1, 0))
dat$lvh <- as.integer(dat$ECG_RV5_SV1_mV >= 4.0)

# ---------- 五项危险因素 ----------
dat$htn      <- as.integer(dat$高血压_综合 == 1)
dat$hypergly <- as.integer(dat$`血糖(GLU)` >= 7.0 | dat$`糖化血红蛋白(HbA1C)` >= 6.5 | dat$自报糖尿病 == 1)
dat$dyslip   <- as.integer(dat$`胆固醇(CHOL)` >= 6.2 | dat$`低密度脂蛋白胆固醇(LDL-CH)` >= 4.1 |
                           dat$`甘油三酯(TG)` >= 2.3 | dat$`高密度脂蛋白胆固醇(HDL-CH)` < 1.0 |
                           dat$自报血脂异常 == 1)
dat$obesity  <- as.integer(dat$中心性肥胖 == 1 | dat$BMI >= 28)
dat$smoke    <- as.integer(dat$当前吸烟 == 1)
dat$cluster  <- dat$htn + dat$hypergly + dat$dyslip + dat$obesity + dat$smoke

# ---------- CKM 分期（AHA 2023，与课题二_表格图片.R 完全一致） ----------
dat$uacr <- dat$尿微量白蛋白 / dat$尿肌酐 * 8840
dat$egfr <- dat$eGFR_池化
dat$adiposity   <- as.integer(dat$BMI >= 24 | dat$中心性肥胖 == 1)
dat$prediabetes <- as.integer((dat$`血糖(GLU)` >= 5.6 & dat$`血糖(GLU)` < 7.0) |
                              (dat$`糖化血红蛋白(HbA1C)` >= 5.7 & dat$`糖化血红蛋白(HbA1C)` < 6.5))
dat$metabolic   <- as.integer(dat$htn == 1 | dat$hypergly == 1 |
                              dat$`甘油三酯(TG)` >= 2.3 | dat$dyslip == 1)
dat$ckd_mod     <- as.integer((dat$egfr >= 30 & dat$egfr < 60) |
                              (dat$uacr >= 30 & dat$uacr < 300))
dat$ckd_high    <- as.integer((dat$egfr >= 15 & dat$egfr < 30) | dat$uacr >= 300)
dat$cvd         <- as.integer(dat$自报心肌梗死 == 1 | dat$自报脑卒中 == 1 |
                              dat$自报冠脉支架 == 1 | dat$自报冠脉搭桥 == 1 |
                              dat$自报房颤 == 1 | dat$自报不稳定心绞痛住院 == 1 |
                              dat$自报外周动脉疾病 == 1)
dat$kidney_fail <- as.integer(dat$egfr < 15)
dat$ckm <- NA_integer_
dat$ckm[dat$cvd == 1 | dat$kidney_fail == 1] <- 4
dat$ckm[is.na(dat$ckm) & dat$ckd_high == 1] <- 3
dat$ckm[is.na(dat$ckm) & (dat$metabolic == 1 | dat$ckd_mod == 1)] <- 2
dat$ckm[is.na(dat$ckm) & (dat$adiposity == 1 | dat$prediabetes == 1)] <- 1
dat$ckm[is.na(dat$ckm)] <- 0
dat$ckm <- factor(dat$ckm, levels = 0:4)

cat("=== CKM 分期分布（核对）===\n"); print(table(dat$ckm))
cat("n =", nrow(dat), "\n\n")

# ============ 敏感性分析 1：排除 Stage 4 后 QTc 趋势 ============
cat("========== [R1-M2] 排除 Stage 4 后 QTc 趋势 ==========\n")
sub <- dat[dat$ckm != "4", ]
sub$ckm_num <- as.numeric(as.character(sub$ckm))
for (s in 0:3) {
  g <- sub[sub$ckm_num == s, ]
  cat(sprintf("  Stage %d: n=%d, QTc事件=%d, %.1f%%\n", s, nrow(g), sum(g$qtc_prol), mean(g$qtc_prol)*100))
}
m_unadj <- glm(qtc_prol ~ ckm_num, data = sub, family = binomial)
m_adj   <- glm(qtc_prol ~ ckm_num + 年龄 + sex, data = sub, family = binomial)
cat(sprintf("\n  排除Stage4后 趋势(未调整): OR=%.2f, P=%.3f\n",
    exp(coef(m_unadj)["ckm_num"]), summary(m_unadj)$coefficients["ckm_num",4]))
cat(sprintf("  排除Stage4后 趋势(调整年龄性别): OR=%.2f, P=%.3f\n",
    exp(coef(m_adj)["ckm_num"]), summary(m_adj)$coefficients["ckm_num",4]))

# 完整 0-4 趋势对照
dat$ckm_num <- as.numeric(as.character(dat$ckm))
m_full_adj <- glm(qtc_prol ~ ckm_num + 年龄 + sex, data = dat, family = binomial)
cat(sprintf("  [对照] 完整0-4趋势(调整年龄性别): OR=%.2f, P=%.3f\n",
    exp(coef(m_full_adj)["ckm_num"]), summary(m_full_adj)$coefficients["ckm_num",4]))

# ============ 敏感性分析 2：CKM vs 聚集度 增量价值 ============
cat("\n========== [R1-M1/R2-M2] CKM vs 聚集度 增量价值 ==========\n")
m_base    <- glm(qtc_prol ~ 年龄 + sex, data = dat, family = binomial)
m_cluster <- glm(qtc_prol ~ 年龄 + sex + cluster, data = dat, family = binomial)
m_ckm     <- glm(qtc_prol ~ 年龄 + sex + ckm_num, data = dat, family = binomial)

lr <- function(fit_full, fit_red) {
  d <- fit_red$deviance - fit_full$deviance
  p <- pchisq(d, 1, lower.tail = FALSE)
  c(LR = d, P = p)
}
lrc <- lr(m_cluster, m_base); lrk <- lr(m_ckm, m_base)
cat(sprintf("  基础模型 -2logLik = %.2f  (AIC=%.1f)\n", m_base$deviance, AIC(m_base)))
cat(sprintf("  +聚集度(连续)  -2logLik = %.2f  LR=%.2f P=%.3f  AIC=%.1f\n",
    m_cluster$deviance, lrc[1], lrc[2], AIC(m_cluster)))
cat(sprintf("  +CKM分期(连续) -2logLik = %.2f  LR=%.2f P=%.3f  AIC=%.1f\n",
    m_ckm$deviance, lrk[1], lrk[2], AIC(m_ckm)))
cat(sprintf("  → CKM AIC vs 聚集度 AIC: %.1f vs %.1f (差值 %.1f, 负值=CKM更优)\n",
    AIC(m_ckm), AIC(m_cluster), AIC(m_ckm)-AIC(m_cluster)))

# AUC（Wilcoxon rank-based）
auc <- function(y, p) {
  n1 <- sum(y == 1); n0 <- sum(y == 0)
  r <- rank(p)
  (sum(r[y == 1]) - n1*(n1+1)/2) / (n1*n0)
}
auc_b <- auc(dat$qtc_prol, predict(m_base, type="response"))
auc_c <- auc(dat$qtc_prol, predict(m_cluster, type="response"))
auc_k <- auc(dat$qtc_prol, predict(m_ckm, type="response"))
cat(sprintf("  AUC 基础=%.3f  +聚集度=%.3f  +CKM=%.3f\n", auc_b, auc_c, auc_k))

# 同时放聚集度+CKM，看 CKM 在控制聚集度后是否仍显著（共线性警告）
m_both <- glm(qtc_prol ~ 年龄 + sex + cluster + ckm_num, data = dat, family = binomial)
cat(sprintf("  同时纳入(聚集度+CKM): cluster P=%.3f, ckm P=%.3f\n",
    summary(m_both)$coefficients["cluster",4],
    summary(m_both)$coefficients["ckm_num",4]))

# ============ 敏感性分析 3：年龄分层 CKM 趋势 ============
cat("\n========== [年龄分层] CKM 分期趋势 ==========\n")
dat$age_grp <- cut(dat$年龄, breaks = c(0, 50, 65, 200), labels = c("<50", "50-64", ">=65"))
for (g in c("<50", "50-64", ">=65")) {
  gg <- dat[dat$age_grp == g, ]
  n <- nrow(gg); k <- sum(gg$qtc_prol)
  if (k < 5 || n < 30) { cat(sprintf("  %s: n=%d, QTc事件=%d（事件过少，跳过）\n", g, n, k)); next }
  m <- glm(qtc_prol ~ ckm_num, data = gg, family = binomial)
  cat(sprintf("  %s: n=%d, QTc事件=%d, CKM趋势 OR=%.2f P=%.3f\n",
      g, n, k, exp(coef(m)["ckm_num"]), summary(m)$coefficients["ckm_num",4]))
}

# ============ 敏感性分析 4：QTc 连续值（非二分）随 CKM 分期的趋势 ============
cat("\n========== [补充] QTc 连续值(ms) 随 CKM 分期 ==========\n")
for (s in 0:4) {
  g <- dat[dat$ckm_num == s, ]
  cat(sprintf("  Stage %d: QTc 中位 %.1f ms (IQR %.1f-%.1f)\n",
      s, median(g$ECG_QTc_ms, na.rm=TRUE),
      quantile(g$ECG_QTc_ms, 0.25, na.rm=TRUE), quantile(g$ECG_QTc_ms, 0.75, na.rm=TRUE)))
}
# Spearman 趋势
sp <- cor.test(dat$ECG_QTc_ms, dat$ckm_num, method="spearman")
cat(sprintf("  Spearman rho=%.3f, P=%.4f\n", sp$estimate, sp$p.value))