# ============================================================
# 课题二 补充分析：阈值二分类 + 复合 ECG 结局
# 1) 聚集度阈值二分类 ≥3 vs <3（QTc 呈阈值效应）
# 2) 复合结局 = QTc延长 | LVH | QRS≥120（提升事件数，回归更稳）
# ============================================================
suppressPackageStartupMessages({
  library(dplyr); library(readr); library(ggplot2)
  library(openxlsx); library(logistf); library(patchwork)
})

DATA <- "D:/龙虾工作空间/2026-08-07-17-00-35/监测数据/选题二_CKM危险因素聚集/课题二_分析数据.csv"
OUT  <- "D:/龙虾工作空间/2026-08-07-17-00-35/监测数据/选题二_CKM危险因素聚集/表格图片"

dat <- read_csv(DATA, show_col_types = FALSE)
dat$sex <- factor(ifelse(dat$性别 == "男", "Male", "Female"), levels = c("Male", "Female"))

# ---------- 结局 ----------
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$qrs_prol <- as.integer(dat$ECG_QRS_ms >= 120)
dat$composite <- as.integer(dat$qtc_prol == 1 | dat$lvh == 1 | dat$qrs_prol == 1)

# ---------- 暴露 ----------
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
dat$cluster_grp <- factor(ifelse(dat$cluster >= 4, "≥4", as.character(dat$cluster)),
                          levels = c("0","1","2","3","≥4"))
dat$cluster_bin <- factor(ifelse(dat$cluster >= 3, "≥3", "<3"), levels = c("<3","≥3"))

# CKM 分期
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("========== 结局事件数 ==========\n")
cat("QTc延长:", sum(dat$qtc_prol, na.rm=TRUE),
    " LVH:", sum(dat$lvh, na.rm=TRUE),
    " QRS≥120:", sum(dat$qrs_prol, na.rm=TRUE),
    " 复合:", sum(dat$composite, na.rm=TRUE),
    sprintf("(%.1f%%)\n", mean(dat$composite, na.rm=TRUE)*100))

# ---------- 关联函数 ----------
firth_or <- function(form, d) {
  fit <- logistf(form, data=d)
  co <- coef(fit); ci <- confint(fit)
  idx <- 2:length(co)
  data.frame(term=names(co)[idx], OR=exp(co[idx]),
             CI=sprintf("%.2f–%.2f", exp(ci[idx,1]), exp(ci[idx,2])),
             P=fit$prob[idx], stringsAsFactors=FALSE)
}
glm_or <- function(form, d) {
  fit <- glm(form, data=d, family=binomial)
  ci <- confint(fit); s <- summary(fit)$coefficients
  data.frame(term=rownames(s)[-1], OR=exp(coef(fit)[-1]),
             CI=sprintf("%.2f–%.2f", exp(ci[-1,1]), exp(ci[-1,2])),
             P=s[-1,4], stringsAsFactors=FALSE)
}
fit_outcome <- function(outcome) {
  d <- dat; d$y <- dat[[outcome]]
  method <- if (sum(d$y, na.rm=TRUE) < 60) firth_or else glm_or
  list(
    bin  = method(y ~ cluster_bin + 年龄 + sex, d),
    grp  = method(y ~ cluster_grp + 年龄 + sex, d),
    ckm  = method(y ~ ckm + 年龄 + sex, d)
  )
}
res <- list(qtc=fit_outcome("qtc_prol"), lvh=fit_outcome("lvh"), composite=fit_outcome("composite"))

cat("\n========== 阈值二分类 ≥3 vs <3 ==========\n")
for (o in names(res)) {
  cat("\n---", o, "---\n")
  b <- res[[o]]$bin; b <- b[grepl("cluster_bin", b$term),]
  print(b)
}

cat("\n========== 复合结局：聚集度分类 + CKM 分期 ==========\n")
print(res$composite$grp)
print(res$composite$ckm)

# ---------- 输出 CSV ----------
write_out <- function(r, name) {
  d <- rbind(
    cbind(r$bin,  exposure="Clustering ≥3 vs <3"),
    cbind(r$grp,  exposure="Clustering (5-level, ref=0)"),
    cbind(r$ckm,  exposure="CKM stage (ref=Stage 0)"))
  write.csv(d, file.path(OUT, name), row.names=FALSE, fileEncoding="UTF-8")
}
write_out(res$qtc,       "Table5_阈值与复合_QTc.csv")
write_out(res$lvh,       "Table6_阈值与复合_LVH.csv")
write_out(res$composite, "Table7_阈值与复合_复合结局.csv")

# ---------- 更新森林图：三结局 × 聚集度5分类 ----------
forest <- rbind(
  cbind(res$qtc$grp,       outcome="QTc prolongation"),
  cbind(res$lvh$grp,       outcome="LVH"),
  cbind(res$composite$grp, outcome="Composite ECG damage"))
forest <- forest %>% filter(grepl("cluster_grp", term)) %>%
  mutate(term = gsub("cluster_grp","", term),
         term = factor(term, levels=rev(c("1","2","3","≥4"))),
         OR = as.numeric(OR),
         lo = as.numeric(sub("–.*","",CI)), hi = as.numeric(sub(".*–","",CI)),
         outcome = factor(outcome, levels=c("QTc prolongation","LVH","Composite ECG damage")))
fig4 <- ggplot(forest, aes(x=OR, y=term)) +
  geom_point(aes(color=outcome), size=2, position=position_dodge(0.6)) +
  geom_linerange(aes(xmin=lo, xmax=hi, color=outcome), position=position_dodge(0.6), size=0.8) +
  geom_vline(xintercept=1, linetype=2, color="grey50") +
  facet_wrap(~outcome) +
  scale_color_manual(values=c("QTc prolongation"="#2E5F8A","LVH"="#C0392B","Composite ECG damage"="#1D9E75")) +
  scale_x_log10() + labs(x="Odds ratio (vs 0 factors, log scale)", y=NULL) +
  theme_bw(base_size=8) + theme(legend.position="none",
    panel.grid.minor=element_blank())
ggsave(file.path(OUT,"Fig4_森林图_三结局.png"), fig4, width=170, height=75, units="mm", dpi=300)
ggsave(file.path(OUT,"Fig4_森林图_三结局.tiff"), fig4, width=170, height=75, units="mm", dpi=300, compression="lzw")

cat("\n========== 完成：补充分析已输出 ==========\n")
print(list.files(OUT, pattern="Table[567]|Fig4_森林图_三结局"))