# ============================================================
# 课题二：CKM 危险因素聚集 + ECG 亚临床表型关联 —— 表格与图片生成
# 目标期刊：BMC Cardiovascular Disorders（图 300 dpi，85/170 mm）
# 结局口径对齐选题一定稿：
#   QTc 延长 = 男≥450 / 女≥460 ms
#   LVH（左室高电压）= Sokolow-Lyon 统一阈值 RV5+SV1 ≥ 4.0 mV（不分性别）
# ============================================================
suppressPackageStartupMessages({
  library(dplyr); library(readr); library(ggplot2)
  library(openxlsx); library(logistf); library(patchwork); library(scales); library(tidyr)
})

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

dat <- read_csv(DATA, show_col_types = FALSE)
dat$sex <- ifelse(dat$性别 == "男", "Male", "Female")
dat$sex <- factor(dat$sex, 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)   # 统一阈值 ≥4.0 mV

# ---------- 五项危险因素 ----------
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)

# 聚集度（0-5），4/5 合并
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"))

# ---------- CKM 分期（AHA 2023 映射） ----------
# UACR (mg/g) ≈ 尿微量白蛋白(mg/L)/尿肌酐(μmol/L) × 8840（单位假设，待核原始单位）
dat$uacr <- dat$尿微量白蛋白 / dat$尿肌酐 * 8840
dat$egfr <- dat$eGFR_池化

dat$adiposity   <- as.integer(dat$BMI >= 24 | dat$中心性肥胖 == 1)                 # Stage1 超重/肥胖
dat$prediabetes <- as.integer((dat$`血糖(GLU)` >= 5.6 & dat$`血糖(GLU)` < 7.0) |
                              (dat$`糖化血红蛋白(HbA1C)` >= 5.7 & dat$`糖化血红蛋白(HbA1C)` < 6.5))  # Stage1 糖前期
dat$metabolic   <- as.integer(dat$htn == 1 | dat$hypergly == 1 |
                              dat$`甘油三酯(TG)` >= 2.3 | dat$dyslip == 1)          # Stage2 代谢
dat$ckd_mod     <- as.integer((dat$egfr >= 30 & dat$egfr < 60) |
                              (dat$uacr >= 30 & dat$uacr < 300))                     # Stage2 中危CKD
dat$ckd_high    <- as.integer((dat$egfr >= 15 & dat$egfr < 30) | dat$uacr >= 300)    # Stage3 高危CKD
dat$cvd         <- as.integer(dat$自报心肌梗死 == 1 | dat$自报脑卒中 == 1 |
                              dat$自报冠脉支架 == 1 | dat$自报冠脉搭桥 == 1 |
                              dat$自报房颤 == 1 | dat$自报不稳定心绞痛住院 == 1 |
                              dat$自报外周动脉疾病 == 1)                             # Stage4 临床CVD
dat$kidney_fail <- as.integer(dat$egfr < 15)                                          # Stage4 肾衰

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)

# 饮酒史：NA 视为未饮酒（对齐选题一）
dat$drink <- as.integer(ifelse(is.na(dat$饮酒史), 0, 1))

cat("========== 关键检查 ==========\n")
cat("QTc延长:", sum(dat$qtc_prol, na.rm=TRUE), "(",
    sprintf("%.1f", mean(dat$qtc_prol, na.rm=TRUE)*100), "%)\n")
cat("LVH(RV5+SV1≥4.0):", sum(dat$lvh, na.rm=TRUE), "(",
    sprintf("%.1f", mean(dat$lvh, na.rm=TRUE)*100), "%)\n")
cat("聚集度分布:\n"); print(table(dat$cluster_grp))
cat("CKM分期分布:\n"); print(table(dat$ckm, useNA="ifany"))
cat("UACR(mg/g) 中位:", sprintf("%.1f", median(dat$uacr, na.rm=TRUE)),
    " 范围[", sprintf("%.1f", min(dat$uacr, na.rm=TRUE)), ",",
    sprintf("%.1f", max(dat$uacr, na.rm=TRUE)), "]\n")
cat("UACR分级: <30:", sum(dat$uacr<30,na.rm=TRUE),
    " 30-300:", sum(dat$uacr>=30&dat$uacr<300,na.rm=TRUE),
    " ≥300:", sum(dat$uacr>=300,na.rm=TRUE), "\n")

# ============================================================
# 统计辅助函数
# ============================================================
# 描述：连续变量中位(IQR)或mean(SD)
med_iqr <- function(x) sprintf("%.1f (%.1f–%.1f)", median(x,na.rm=TRUE),
                               quantile(x,0.25,na.rm=TRUE), quantile(x,0.75,na.rm=TRUE))
n_pct <- function(x) { t <- table(x); paste0(t, " (", sprintf("%.1f", t/sum(t)*100), ")") }

# Firth 惩罚回归提取 OR/CI/P（用于 LVH 稀疏事件）
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
  )
}
# 普通 logistic 提取 OR/CI/P
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
  )
}

# ============================================================
# Table 1：按 CKM 分期分层的基线特征
# ============================================================
vars_cont <- c("年龄","BMI","SBP均值_左23","DBP均值_左23","血糖(GLU)","糖化血红蛋白(HbA1C)",
               "低密度脂蛋白胆固醇(LDL-CH)","甘油三酯(TG)","高密度脂蛋白胆固醇(HDL-CH)",
               "eGFR_池化","uacr","ECG_QTc_ms","ECG_RV5_SV1_mV")
t1 <- lapply(split(dat, dat$ckm), function(g) {
  c(
    n = nrow(g),
    Male = sum(g$sex=="Male",na.rm=TRUE),
    age = med_iqr(g$年龄),
    bmi = med_iqr(g$BMI),
    sbp = med_iqr(g$SBP均值_左23),
    glu = med_iqr(g$`血糖(GLU)`),
    hba = med_iqr(g$`糖化血红蛋白(HbA1C)`),
    ldl = med_iqr(g$`低密度脂蛋白胆固醇(LDL-CH)`),
    tg  = med_iqr(g$`甘油三酯(TG)`),
    hdl = med_iqr(g$`高密度脂蛋白胆固醇(HDL-CH)`),
    egfr = med_iqr(g$eGFR_池化),
    uacr = med_iqr(g$uacr),
    qtc = med_iqr(g$ECG_QTc_ms),
    rv5sv1 = med_iqr(g$ECG_RV5_SV1_mV),
    qtc_prol = sum(g$qtc_prol,na.rm=TRUE),
    lvh = sum(g$lvh,na.rm=TRUE)
  )
})
T1 <- as.data.frame(t1, check.names=FALSE)
T1 <- cbind(`Variable`=c("N","Male, n (%)","Age, y","BMI, kg/m²","SBP, mmHg","Glucose, mmol/L",
                         "HbA1c, %","LDL-C, mmol/L","TG, mmol/L","HDL-C, mmol/L",
                         "eGFR, mL/min/1.73m²","UACR, mg/g","QTc, ms","RV5+SV1, mV",
                         "QTc prolongation, n (%)","LVH, n (%)"), T1)
write.csv(T1, file.path(OUT,"Table1_CKM基线.csv"), row.names=FALSE, fileEncoding="UTF-8")

# ============================================================
# Table 2：聚集度分布（总体 + 性别）
# ============================================================
T2 <- dat %>% group_by(cluster_grp) %>% summarise(
  Total = n(),
  Male = sum(sex=="Male",na.rm=TRUE),
  Female = sum(sex=="Female",na.rm=TRUE)
) %>% mutate(Total_pct = sprintf("%.1f", Total/sum(Total)*100))
write.csv(T2, file.path(OUT,"Table2_聚集度分布.csv"), row.names=FALSE, fileEncoding="UTF-8")

# ============================================================
# Table 3 & 4：聚集度 / CKM 分期 → QTc、LVH 关联
# ============================================================
# 趋势检验用连续编码
run_assoc <- function(outcome, label) {
  d <- dat; d$y <- dat[[outcome]]
  # 聚集度（分类 + 趋势）
  m_cat <- if (outcome=="lvh") {
    firth_or(y ~ cluster_grp + 年龄 + sex, d)
  } else {
    glm_or(y ~ cluster_grp + 年龄 + sex, d)
  }
  m_trend <- if (outcome=="lvh") {
    firth_or(y ~ cluster + 年龄 + sex, d)
  } else {
    glm_or(y ~ cluster + 年龄 + sex, d)
  }
  # CKM 分期（分类 + 趋势）
  c_cat <- if (outcome=="lvh") {
    firth_or(y ~ ckm + 年龄 + sex, d)
  } else {
    glm_or(y ~ ckm + 年龄 + sex, d)
  }
  c_trend <- if (outcome=="lvh") {
    firth_or(y ~ as.numeric(ckm) + 年龄 + sex, d)
  } else {
    glm_or(y ~ as.numeric(ckm) + 年龄 + sex, d)
  }
  list(cluster_cat=m_cat, cluster_trend=m_trend, ckm_cat=c_cat, ckm_trend=c_trend)
}
res_qtc <- run_assoc("qtc_prol","QTc")
res_lvh <- run_assoc("lvh","LVH")

# 汇总关联表
bind_assoc <- function(res, lab) {
  a <- res$cluster_cat; a$exposure <- "Clustering (ref=0)"
  b <- res$ckm_cat; b$exposure <- "CKM stage (ref=Stage 0)"
  tt <- res$cluster_trend; tt$exposure <- "Clustering (per +1 factor)"
  ct <- res$ckm_trend; ct$exposure <- "CKM stage (per +1 stage)"
  rbind(a,b,tt,ct)
}
T3 <- bind_assoc(res_qtc,"QTc"); T3$Outcome <- "QTc prolongation"
T4 <- bind_assoc(res_lvh,"LVH"); T4$Outcome <- "LVH"
write.csv(T3, file.path(OUT,"Table3_聚集度_CKM_与QTc关联.csv"), row.names=FALSE, fileEncoding="UTF-8")
write.csv(T4, file.path(OUT,"Table4_聚集度_CKM_与LVH关联.csv"), row.names=FALSE, fileEncoding="UTF-8")

# 打印关联结果供核对
cat("\n========== QTc 关联（聚集度 分类）==========\n"); print(res_qtc$cluster_cat)
cat("\n========== QTc 关联（CKM 分类）==========\n"); print(res_qtc$ckm_cat)
cat("\n========== LVH 关联（Firth，聚集度 分类）==========\n"); print(res_lvh$cluster_cat)
cat("\n========== LVH 关联（Firth，CKM 分类）==========\n"); print(res_lvh$ckm_cat)

# ============================================================
# 图片（BMC 规范：300 dpi，单栏 85mm / 双栏 170mm）
# ============================================================
theme_set(theme_bw(base_size=8) + theme(
  panel.grid.minor=element_blank(),
  plot.title=element_text(size=9, face="bold"),
  axis.title=element_text(size=8), axis.text=element_text(size=7)))

# Fig 1：聚集度分布柱状图
fig1 <- ggplot(dat, aes(x=cluster_grp)) + geom_bar(fill="#2E5F8A", width=0.7) +
  geom_text(stat="count", aes(label=after_stat(count)), vjust=-0.4, size=2.5) +
  labs(x="Number of clustered risk factors", y="Participants (n)") +
  scale_y_continuous(expand=expansion(mult=c(0,0.08)))
ggsave(file.path(OUT,"Fig1_聚集度分布.png"), fig1, width=85, height=70, units="mm", dpi=300)
ggsave(file.path(OUT,"Fig1_聚集度分布.tiff"), fig1, width=85, height=70, units="mm", dpi=300, compression="lzw")

# Fig 2：CKM 分期构成
fig2 <- ggplot(dat, aes(x=ckm)) + geom_bar(fill="#C0392B", width=0.7) +
  geom_text(stat="count", aes(label=after_stat(count)), vjust=-0.4, size=2.5) +
  labs(x="CKM stage (AHA 2023)", y="Participants (n)") +
  scale_y_continuous(expand=expansion(mult=c(0,0.08)))
ggsave(file.path(OUT,"Fig2_CKM分期.png"), fig2, width=85, height=70, units="mm", dpi=300)
ggsave(file.path(OUT,"Fig2_CKM分期.tiff"), fig2, width=85, height=70, units="mm", dpi=300, compression="lzw")

# Fig 3：QTc/LVH 患病率随聚集度趋势
prev <- dat %>% group_by(cluster_grp) %>% summarise(
  QTc = mean(qtc_prol,na.rm=TRUE)*100, LVH = mean(lvh,na.rm=TRUE)*100) %>%
  tidyr::pivot_longer(-cluster_grp, names_to="outcome", values_to="prev")
fig3 <- ggplot(prev, aes(x=cluster_grp, y=prev, fill=outcome)) +
  geom_col(position=position_dodge(0.8), width=0.7) +
  geom_text(aes(label=sprintf("%.1f",prev)), position=position_dodge(0.8), vjust=-0.4, size=2.3) +
  scale_fill_manual(values=c("QTc"="#2E5F8A","LVH"="#C0392B")) +
  labs(x="Number of clustered risk factors", y="Prevalence (%)", fill=NULL) +
  theme(legend.position="top") + scale_y_continuous(expand=expansion(mult=c(0,0.1)))
ggsave(file.path(OUT,"Fig3_患病率趋势.png"), fig3, width=170, height=80, units="mm", dpi=300)
ggsave(file.path(OUT,"Fig3_患病率趋势.tiff"), fig3, width=170, height=80, units="mm", dpi=300, compression="lzw")

# Fig 4：森林图（聚集度各水平 vs 0 的 OR，QTc & LVH）
forest_df <- rbind(
  cbind(res_qtc$cluster_cat, outcome="QTc prolongation"),
  cbind(res_lvh$cluster_cat, outcome="LVH"))
forest_df <- forest_df %>% filter(!grepl("年龄|sex|Sex", term)) %>%
  mutate(term = gsub("cluster_grp","",term),
         term = factor(term, levels=rev(c("1","2","3","≥4"))),
         OR = as.numeric(OR))
fig4 <- ggplot(forest_df, aes(x=OR, y=term, xmin=1, xmax=1)) +
  geom_point(aes(color=outcome), size=2) +
  geom_errorbarh(aes(xmin=as.numeric(sub("–.*","",CI)), xmax=as.numeric(sub(".*–","",CI)),
                     color=outcome), height=0.2) +
  geom_vline(xintercept=1, linetype=2, color="grey50") +
  facet_wrap(~outcome) +
  scale_color_manual(values=c("QTc prolongation"="#2E5F8A","LVH"="#C0392B")) +
  scale_x_log10() + labs(x="Odds ratio (vs 0 factors, log scale)", y=NULL) +
  theme(legend.position="none")
ggsave(file.path(OUT,"Fig4_森林图.png"), fig4, width=170, height=70, units="mm", dpi=300)
ggsave(file.path(OUT,"Fig4_森林图.tiff"), fig4, width=170, height=70, units="mm", dpi=300, compression="lzw")

cat("\n========== 完成 ==========\n")
cat("表格与图片已输出至:", OUT, "\n")
cat("文件清单:\n"); print(list.files(OUT))