# ============================================================
# 课题二 核心结果表重排（重新定位版）
# 定位：描述性贡献（CKM分期/聚集度分布）+ QTc延长为主结局 + 聚集度5分类(不合并) + CKM分期趋势
# ============================================================
suppressPackageStartupMessages({
  library(dplyr); library(readr); library(openxlsx); library(logistf)
})

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$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$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_num <- dat$ckm
dat$ckm <- factor(dat$ckm, levels = 0:4)

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))
trend_p_cont <- function(x) {   # 有序分期 vs 连续变量：Spearman 秩相关趋势
  cc <- complete.cases(x, dat$ckm_num)
  if (sum(cc) < 3) return("")
  p <- suppressWarnings(cor.test(x[cc], dat$ckm_num[cc], method="spearman")$p.value)
  fmt_p(p)
}
trend_p_cat <- function(x) {    # 有序分期 vs 二分类：Cochran-Armitage
  tt <- table(x, dat$ckm_num)
  if (nrow(tt) < 2) return("")
  fmt_p(prop.trend.test(tt[2,], colSums(tt))$p.value)
}
fmt_p <- function(p) {
  p <- suppressWarnings(as.numeric(p))
  ifelse(is.na(p), "", ifelse(p < 0.001, "<0.001", sprintf("%.3f", p)))
}
fmt_or <- function(x) sprintf("%.2f", as.numeric(x))

# ============================================================
# Table 1：按 CKM 分期分层的基线特征（median(IQR) + 趋势P）
# ============================================================
stages <- as.character(0:4)
row_build <- list()
for (v in c("age","bmi","sbp","dbp","glu","hba","ldl","tg","hdl","egfr","uacr","qtc","rv5sv1")) {
  d <- dat
  d$x <- switch(v,
    age=d$年龄, bmi=d$BMI, sbp=d$SBP均值_左23, dbp=d$DBP均值_左23,
    glu=d$`血糖(GLU)`, hba=d$`糖化血红蛋白(HbA1C)`, ldl=d$`低密度脂蛋白胆固醇(LDL-CH)`,
    tg=d$`甘油三酯(TG)`, hdl=d$`高密度脂蛋白胆固醇(HDL-CH)`, egfr=d$eGFR_池化,
    uacr=d$uacr, qtc=d$ECG_QTc_ms, rv5sv1=d$ECG_RV5_SV1_mV)
  row_build[[v]] <- c(sapply(stages, function(s) med_iqr(d$x[d$ckm==s])),
                      med_iqr(d$x), trend_p_cont(d$x))}
cat_row <- function(v, x) {
  tt <- table(x, dat$ckm_num)
  c(sprintf("%d (%s)", tt[2,], paste0(sprintf("%.1f", tt[2,]/colSums(tt)*100), "%")),
    sprintf("%d (%.1f%%)", sum(tt[2,]), sum(tt[2,])/sum(tt)*100),
    trend_p_cat(x))
}
T1 <- rbind(
  `N` = c(sapply(stages, function(s) sum(dat$ckm==s)), nrow(dat), NA),
  `Male, n (%)` = cat_row("male", dat$sex=="Male"),
  `Age, years` = row_build$age,
  `BMI, kg/m²` = row_build$bmi,
  `SBP, mmHg` = row_build$sbp,
  `DBP, mmHg` = row_build$dbp,
  `Fasting glucose, mmol/L` = row_build$glu,
  `HbA1c, %` = row_build$hba,
  `LDL-C, mmol/L` = row_build$ldl,
  `Triglycerides, mmol/L` = row_build$tg,
  `HDL-C, mmol/L` = row_build$hdl,
  `eGFR, mL/min/1.73m²` = row_build$egfr,
  `UACR, mg/g` = row_build$uacr,
  `QTc, ms` = row_build$qtc,
  `RV5+SV1, mV` = row_build$rv5sv1,
  `QTc prolongation, n (%)` = cat_row("qtc", dat$qtc_prol),
  `LVH, n (%)` = cat_row("lvh", dat$lvh),
  `QRS ≥120 ms, n (%)` = cat_row("qrs", dat$qrs_prol)
)
colnames(T1) <- c(paste0("Stage ", stages), "Total", "P for trend")
T1 <- data.frame(Variable=rownames(T1), T1, check.names=FALSE, row.names=NULL)

# ============================================================
# Table 2：危险因素单项患病率 + 聚集度分布 + CKM 分期分布
# ============================================================
rf <- data.frame(
  Risk_factor = c("Hypertension","Hyperglycemia","Dyslipidemia","Obesity","Smoking"),
  n = c(sum(dat$htn,na.rm=TRUE),sum(dat$hypergly,na.rm=TRUE),sum(dat$dyslip,na.rm=TRUE),
        sum(dat$obesity,na.rm=TRUE),sum(dat$smoke,na.rm=TRUE)),
  stringsAsFactors = FALSE)
rf$pct <- sprintf("%.1f", rf$n/nrow(dat)*100)
clust_tab <- as.data.frame(table(dat$cluster_grp), stringsAsFactors=FALSE)
names(clust_tab) <- c("Clustering","n")
clust_tab$pct <- sprintf("%.1f", clust_tab$n/nrow(dat)*100)
ckm_tab <- as.data.frame(table(dat$ckm), stringsAsFactors=FALSE)
names(ckm_tab) <- c("CKM_stage","n")
ckm_tab$pct <- sprintf("%.1f", ckm_tab$n/nrow(dat)*100)

# ============================================================
# Table 3 & 4：核心关联（QTc 主 / LVH 次），聚集度5分类 + CKM分期 + 趋势
# ============================================================
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)
}
assoc_table <- function(outcome, method) {
  d <- dat; d$y <- dat[[outcome]]
  m <- if (method=="firth") firth_or else glm_or
  g <- m(y ~ cluster_grp + 年龄 + sex, d)
  c <- m(y ~ ckm + 年龄 + sex, d)
  gt <- m(y ~ cluster + 年龄 + sex, d)
  ct <- m(y ~ as.numeric(ckm) + 年龄 + sex, d)
  # 提取暴露相关行，重排成清晰格式
  pick <- function(r, pat) r[grepl(pat, r$term), c("term","OR","CI","P")]
  gc <- pick(g, "cluster_grp"); gc$term <- gsub("cluster_grp","", gc$term)
  cc <- pick(c, "ckm");        cc$term <- gsub("ckm","Stage ", cc$term)
  gt <- pick(gt, "^cluster$"); gt$term <- "Per +1 factor"
  ct <- pick(ct, "as.numeric"); ct$term <- "Per +1 stage"
  gc$OR <- fmt_or(gc$OR); gc$P <- fmt_p(gc$P)
  cc$OR <- fmt_or(cc$OR); cc$P <- fmt_p(cc$P)
  gt$OR <- fmt_or(gt$OR); gt$P <- fmt_p(gt$P)
  ct$OR <- fmt_or(ct$OR); ct$P <- fmt_p(ct$P)
  rbind(
    data.frame(Exposure="Clustering (ref = 0 factors)", Level=gc$term, OR=gc$OR, `95% CI`=gc$CI, P=gc$P, check.names=FALSE),
    data.frame(Exposure="", Level=gt$term, OR=gt$OR, `95% CI`=gt$CI, P=gt$P, check.names=FALSE),
    data.frame(Exposure="CKM stage (ref = Stage 0)", Level=cc$term, OR=cc$OR, `95% CI`=cc$CI, P=cc$P, check.names=FALSE),
    data.frame(Exposure="", Level=ct$term, OR=ct$OR, `95% CI`=ct$CI, P=ct$P, check.names=FALSE)
  )
}
T3 <- assoc_table("qtc_prol", "glm")     # QTc 主结局，普通 Logistic
T4 <- assoc_table("lvh", "firth")        # LVH 次结局，Firth

# ============================================================
# 输出 CSV + 汇总 xlsx
# ============================================================
write.csv(T1, file.path(OUT,"核心表1_CKM基线.csv"), row.names=FALSE, fileEncoding="UTF-8")
write.csv(rf, file.path(OUT,"核心表2a_危险因素.csv"), row.names=FALSE, fileEncoding="UTF-8")
write.csv(clust_tab, file.path(OUT,"核心表2b_聚集度分布.csv"), row.names=FALSE, fileEncoding="UTF-8")
write.csv(ckm_tab, file.path(OUT,"核心表2c_CKM分期分布.csv"), row.names=FALSE, fileEncoding="UTF-8")
write.csv(T3, file.path(OUT,"核心表3_QTc关联.csv"), row.names=FALSE, fileEncoding="UTF-8")
write.csv(T4, file.path(OUT,"核心表4_LVH关联_探索.csv"), row.names=FALSE, fileEncoding="UTF-8")

wb <- createWorkbook()
addWorksheet(wb, "Table1_CKM基线");  writeData(wb, "Table1_CKM基线", T1)
addWorksheet(wb, "Table2_危险因素"); writeData(wb, "Table2_危险因素", rf)
addWorksheet(wb, "Table2_聚集度");   writeData(wb, "Table2_聚集度", clust_tab)
addWorksheet(wb, "Table2_CKM分期");  writeData(wb, "Table2_CKM分期", ckm_tab)
addWorksheet(wb, "Table3_QTc关联");  writeData(wb, "Table3_QTc关联", T3)
addWorksheet(wb, "Table4_LVH关联");  writeData(wb, "Table4_LVH关联", T4)
saveWorkbook(wb, file.path(OUT, "核心结果表_重排版.xlsx"), overwrite=TRUE)

cat("========== 核心表1：CKM基线（median IQR + 趋势P）==========\n")
print(T1, row.names=FALSE, width=200)
cat("\n========== 核心表3：QTc关联（主结局）==========\n")
print(T3, row.names=FALSE)
cat("\n========== 核心表4：LVH关联（探索）==========\n")
print(T4, row.names=FALSE)
cat("\n已完成，输出至:", OUT, "\n")