## =========================
## 0. 环境准备
## =========================

rm(list = ls())
options(stringsAsFactors = FALSE)
set.seed(2025)

work_dir <- "C:/Users/lenovo/Documents/肌间静脉血栓"
setwd(work_dir)

out_dir <- file.path(work_dir, "analysis_outputs")
if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)

install_if_missing <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
  }
}

pkgs <- c(
  "readxl", "openxlsx", "dplyr", "tidyr", "moments",
  "glmnet", "rms", "pROC", "ResourceSelection"
)

invisible(lapply(pkgs, install_if_missing))
invisible(lapply(pkgs, library, character.only = TRUE))

parse_num <- function(x) {
  as.numeric(gsub("[^0-9eE+\\.-]", "", as.character(x)))
}

fmt_p <- function(p) {
  ifelse(is.na(p), NA,
         ifelse(p < 0.001, "<0.001", sprintf("%.3f", p)))
}

safe_logit <- function(p) {
  p <- pmin(pmax(p, 1e-6), 1 - 1e-6)
  log(p / (1 - p))
}

get_auc <- function(y, p) {
  as.numeric(pROC::auc(pROC::roc(y, p, quiet = TRUE, direction = "<")))
}


## =========================
## 1. 读取并整理数据
## =========================

df <- readxl::read_excel(file.path(work_dir, "cleaned_data.xlsx"))
df <- as.data.frame(df)

## 将 TG 名称统一为 TG；原始数据列名为 Triglycerides
if ("Triglycerides" %in% names(df) && !"TG" %in% names(df)) {
  df$TG <- df$Triglycerides
}

## 结局变量：MCVT，有 = 1，无 = 0
df$MCVT <- trimws(as.character(df$MCVT))
df$MCVT_bin <- ifelse(df$MCVT == "有", 1, 0)
df$MCVT_group <- factor(df$MCVT_bin, levels = c(0, 1),
                        labels = c("Non-MCVT", "MCVT"))

continuous_vars <- c(
  "Age", "BMI", "Time", "CCI", "WBC", "N", "L", "Singlecore",
  "HB", "HCT", "PLT", "NLR", "PLR", "TP", "Alb", "Tbil",
  "Crea", "Urea", "TG", "TC", "EPI", "P", "INR", "Fbg",
  "APTT", "TT", "PT", "DDI"
)

categorical_vars <- c(
  "Gender", "Hypertension", "Diabate", "Coronary", "Fracture_type"
)

continuous_vars <- intersect(continuous_vars, names(df))
categorical_vars <- intersect(categorical_vars, names(df))

df[continuous_vars] <- lapply(df[continuous_vars], parse_num)

for (v in categorical_vars) {
  df[[v]] <- as.factor(df[[v]])
}


## =========================
## 2. 所有变量单因素分析
## =========================

check_normality <- function(x) {
  x <- x[!is.na(x)]
  if (length(x) < 3) return(FALSE)
  sk <- moments::skewness(x)
  ku <- moments::kurtosis(x)
  if (is.na(sk) || is.na(ku)) return(FALSE)
  abs(sk) <= 2 && ku <= 7
}

mean_sd <- function(x) {
  sprintf("%.2f ± %.2f", mean(x, na.rm = TRUE), sd(x, na.rm = TRUE))
}

median_iqr <- function(x) {
  q <- quantile(x, probs = c(0.25, 0.50, 0.75), na.rm = TRUE)
  sprintf("%.2f (%.2f-%.2f)", q[2], q[1], q[3])
}

univ_list <- list()
univ_p <- data.frame(Variable = character(), P_value = numeric())

for (v in continuous_vars) {
  x0 <- df[[v]][df$MCVT_bin == 0]
  x1 <- df[[v]][df$MCVT_bin == 1]
  x_all <- df[[v]]

  normal_flag <- check_normality(x_all)

  if (normal_flag) {
    overall <- mean_sd(x_all)
    non_mcv <- mean_sd(x0)
    mcv <- mean_sd(x1)
    p_val <- tryCatch(t.test(x0, x1, var.equal = FALSE)$p.value,
                      error = function(e) NA_real_)
    test_name <- "Welch t test"
  } else {
    overall <- median_iqr(x_all)
    non_mcv <- median_iqr(x0)
    mcv <- median_iqr(x1)
    p_val <- tryCatch(wilcox.test(x0, x1, exact = FALSE)$p.value,
                      error = function(e) NA_real_)
    test_name <- "Wilcoxon rank-sum test"
  }

  univ_list[[v]] <- data.frame(
    Variable = v,
    Level = "",
    Overall = overall,
    Non_MCVT = non_mcv,
    MCVT = mcv,
    Test = test_name,
    P_value = p_val,
    P_display = fmt_p(p_val)
  )

  univ_p <- rbind(univ_p, data.frame(Variable = v, P_value = p_val))
}

for (v in categorical_vars) {
  tab <- table(df[[v]], df$MCVT_group, useNA = "no")
  tab <- tab[rowSums(tab) > 0, , drop = FALSE]

  if (nrow(tab) == 0) next

  expected_ok <- FALSE
  chi_p <- NA_real_
  if (nrow(tab) >= 2 && ncol(tab) == 2) {
    chi_obj <- suppressWarnings(chisq.test(tab, correct = FALSE))
    expected_ok <- all(chi_obj$expected >= 5)
    chi_p <- chi_obj$p.value
  }

  if (expected_ok) {
    p_val <- chi_p
    test_name <- "Chi-square test"
  } else {
    p_val <- tryCatch(
      fisher.test(tab)$p.value,
      error = function(e) fisher.test(tab, simulate.p.value = TRUE, B = 100000)$p.value
    )
    if (v == "Fracture_type") {
      test_name <- "Fisher-Freeman-Halton exact test"
    } else {
      test_name <- "Fisher exact test"
    }
  }

  total_n <- sum(tab)
  non_n <- sum(tab[, "Non-MCVT"])
  mcv_n <- sum(tab[, "MCVT"])

  for (i in seq_len(nrow(tab))) {
    level_name <- rownames(tab)[i]
    overall <- sprintf("%d (%.1f%%)", sum(tab[i, ]), 100 * sum(tab[i, ]) / total_n)
    non_mcv <- sprintf("%d (%.1f%%)", tab[i, "Non-MCVT"], 100 * tab[i, "Non-MCVT"] / non_n)
    mcv <- sprintf("%d (%.1f%%)", tab[i, "MCVT"], 100 * tab[i, "MCVT"] / mcv_n)

    univ_list[[paste(v, level_name, sep = "_")]] <- data.frame(
      Variable = ifelse(i == 1, v, ""),
      Level = level_name,
      Overall = overall,
      Non_MCVT = non_mcv,
      MCVT = mcv,
      Test = ifelse(i == 1, test_name, ""),
      P_value = ifelse(i == 1, p_val, NA_real_),
      P_display = ifelse(i == 1, fmt_p(p_val), "")
    )
  }

  univ_p <- rbind(univ_p, data.frame(Variable = v, P_value = p_val))
}

univariate_table <- dplyr::bind_rows(univ_list)

write.csv(univariate_table,
          file.path(out_dir, "univariate_all_variables_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")

write.csv(univ_p,
          file.path(out_dir, "univariate_p_values_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")


## =========================
## 3. LASSO 变量筛选
## =========================

## 筛选规则：
## 单因素分析 P < 0.01 的变量进入 LASSO；
## APTT 基于临床意义和原始建模策略保留。

candidate_vars <- univ_p$Variable[!is.na(univ_p$P_value) & univ_p$P_value < 0.01]
candidate_vars <- unique(c(candidate_vars, "APTT"))
candidate_vars <- intersect(candidate_vars, c(continuous_vars, categorical_vars))

## 如需严格复现返修稿中的候选变量，可使用下面这一行固定候选变量：
candidate_vars <- c("Time", "N", "L", "PLT", "NLR", "PLR", "TG", "TC", "APTT", "Fracture_type")
candidate_vars <- intersect(candidate_vars, names(df))

lasso_df <- df[, c("MCVT_bin", candidate_vars)]
lasso_df <- na.omit(lasso_df)

lasso_formula <- as.formula(
  paste("MCVT_bin ~", paste(candidate_vars, collapse = " + "))
)

x_lasso <- model.matrix(lasso_formula, data = lasso_df)[, -1, drop = FALSE]
y_lasso <- lasso_df$MCVT_bin

set.seed(2025)
cvfit <- glmnet::cv.glmnet(
  x = x_lasso,
  y = y_lasso,
  family = "binomial",
  alpha = 1,
  nfolds = 10,
  standardize = TRUE,
  type.measure = "deviance"
)

lambda_min <- cvfit$lambda.min
lambda_1se <- cvfit$lambda.1se

coef_1se <- as.matrix(coef(cvfit, s = "lambda.1se"))
lasso_coef_table <- data.frame(
  Feature = rownames(coef_1se),
  Coefficient_lambda_1se = as.numeric(coef_1se),
  Selected = as.numeric(coef_1se) != 0
)

selected_lasso <- lasso_coef_table$Feature[
  lasso_coef_table$Selected & lasso_coef_table$Feature != "(Intercept)"
]

write.csv(lasso_coef_table,
          file.path(out_dir, "lasso_coefficients_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")


## =========================
## 4. RCS 分析
## =========================

## TG 使用 3 个 knots，位于第 10、50、90 百分位
tg_knots <- as.numeric(quantile(df$TG, probs = c(0.10, 0.50, 0.90), na.rm = TRUE))

## 若希望完全采用返修稿报告的固定 knot 位置，可取消下一行注释：
## tg_knots <- c(0.62, 1.03, 1.64)

rcs_mat <- rms::rcspline.eval(df$TG, knots = tg_knots, inclx = TRUE)
rcs_mat <- as.data.frame(rcs_mat)
names(rcs_mat) <- c("TG_linear", "TG_rcs1")

df_model <- cbind(df, rcs_mat)

fit_no_tg <- glm(
  MCVT_bin ~ Time + NLR + PLR + TC + APTT,
  data = df_model,
  family = binomial()
)

fit_tg_linear <- glm(
  MCVT_bin ~ Time + NLR + PLR + TG_linear + TC + APTT,
  data = df_model,
  family = binomial()
)

fit_final <- glm(
  MCVT_bin ~ Time + NLR + PLR + TG_linear + TG_rcs1 + TC + APTT,
  data = df_model,
  family = binomial()
)

lr_overall <- 2 * (as.numeric(logLik(fit_final)) - as.numeric(logLik(fit_no_tg)))
df_overall <- attr(logLik(fit_final), "df") - attr(logLik(fit_no_tg), "df")
p_overall <- pchisq(lr_overall, df = df_overall, lower.tail = FALSE)

lr_nonlinear <- 2 * (as.numeric(logLik(fit_final)) - as.numeric(logLik(fit_tg_linear)))
df_nonlinear <- attr(logLik(fit_final), "df") - attr(logLik(fit_tg_linear), "df")
p_nonlinear <- pchisq(lr_nonlinear, df = df_nonlinear, lower.tail = FALSE)

rcs_result <- data.frame(
  Variable = "TG",
  Knots = paste(round(tg_knots, 3), collapse = ", "),
  Overall_LR_chisq = lr_overall,
  Overall_df = df_overall,
  Overall_P = p_overall,
  Nonlinear_LR_chisq = lr_nonlinear,
  Nonlinear_df = df_nonlinear,
  Nonlinear_P = p_nonlinear
)

write.csv(rcs_result,
          file.path(out_dir, "rcs_nonlinearity_results_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")


## =========================
## 5. 最终模型公式和回归系数
## =========================

coef_final <- summary(fit_final)$coefficients

coef_table <- data.frame(
  Term = rownames(coef_final),
  Coefficient_beta = coef_final[, "Estimate"],
  SE = coef_final[, "Std. Error"],
  Wald_z = coef_final[, "z value"],
  P_value = coef_final[, "Pr(>|z|)"],
  OR = exp(coef_final[, "Estimate"]),
  OR_95CI_lower = exp(coef_final[, "Estimate"] - 1.96 * coef_final[, "Std. Error"]),
  OR_95CI_upper = exp(coef_final[, "Estimate"] + 1.96 * coef_final[, "Std. Error"])
)

write.csv(coef_table,
          file.path(out_dir, "full_model_coefficients_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")

## 模型应用：
## LP = intercept + beta1*Time + beta2*NLR + beta3*PLR
##      + beta4*TG_linear + beta5*TG_rcs1 + beta6*TC + beta7*APTT
## Predicted risk = 1 / (1 + exp(-LP))

pred_prob <- as.numeric(predict(fit_final, type = "response"))
lp <- as.numeric(predict(fit_final, type = "link"))


## =========================
## 6. EPV 和 uniform shrinkage factor
## =========================

event_n <- sum(df_model$MCVT_bin == 1, na.rm = TRUE)
non_event_n <- sum(df_model$MCVT_bin == 0, na.rm = TRUE)

fit_null <- glm(MCVT_bin ~ 1, data = df_model, family = binomial())

model_df <- attr(logLik(fit_final), "df") - attr(logLik(fit_null), "df")
lr_chisq <- 2 * (as.numeric(logLik(fit_final)) - as.numeric(logLik(fit_null)))

parameter_based_epv <- event_n / model_df
uniform_shrinkage_factor <- (lr_chisq - model_df) / lr_chisq

epv_table <- data.frame(
  Events = event_n,
  Non_events = non_event_n,
  Predictor_parameters_excluding_intercept = model_df,
  TG_spline_parameters = 2,
  Parameter_based_EPV = parameter_based_epv,
  LR_chisq = lr_chisq,
  Uniform_shrinkage_factor = uniform_shrinkage_factor
)

write.csv(epv_table,
          file.path(out_dir, "epv_shrinkage_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")


## =========================
## 7. AUC、Brier score、C-index 及 bootstrap CI
## =========================

auc_app <- get_auc(df_model$MCVT_bin, pred_prob)
auc_ci <- as.numeric(pROC::ci.auc(
  pROC::roc(df_model$MCVT_bin, pred_prob, quiet = TRUE, direction = "<")
))

brier_app <- mean((pred_prob - df_model$MCVT_bin)^2)
c_index_app <- auc_app

bootstrap_model_metrics <- function(data, indices) {
  d_boot <- data[indices, ]

  fit_b <- tryCatch(
    glm(MCVT_bin ~ Time + NLR + PLR + TG_linear + TG_rcs1 + TC + APTT,
        data = d_boot, family = binomial()),
    error = function(e) NULL
  )

  if (is.null(fit_b)) {
    return(c(
      auc_boot = NA, auc_test = NA,
      brier_boot = NA, brier_test = NA,
      c_boot = NA, c_test = NA
    ))
  }

  p_boot <- as.numeric(predict(fit_b, newdata = d_boot, type = "response"))
  p_test <- as.numeric(predict(fit_b, newdata = data, type = "response"))

  auc_boot <- tryCatch(get_auc(d_boot$MCVT_bin, p_boot), error = function(e) NA_real_)
  auc_test <- tryCatch(get_auc(data$MCVT_bin, p_test), error = function(e) NA_real_)

  brier_boot <- mean((p_boot - d_boot$MCVT_bin)^2)
  brier_test <- mean((p_test - data$MCVT_bin)^2)

  c(auc_boot = auc_boot, auc_test = auc_test,
    brier_boot = brier_boot, brier_test = brier_test,
    c_boot = auc_boot, c_test = auc_test)
}

B <- 1000
set.seed(2025)
boot_indices <- replicate(B, sample(seq_len(nrow(df_model)), replace = TRUE))

boot_metrics <- t(apply(boot_indices, 2, function(idx) {
  bootstrap_model_metrics(df_model, idx)
}))

auc_optimism <- mean(boot_metrics[, "auc_boot"] - boot_metrics[, "auc_test"], na.rm = TRUE)
brier_optimism <- mean(boot_metrics[, "brier_boot"] - boot_metrics[, "brier_test"], na.rm = TRUE)
c_optimism <- mean(boot_metrics[, "c_boot"] - boot_metrics[, "c_test"], na.rm = TRUE)

auc_corrected <- auc_app - auc_optimism
brier_corrected <- brier_app - brier_optimism
c_index_corrected <- c_index_app - c_optimism

brier_boot_values <- boot_metrics[, "brier_test"]
c_boot_values <- boot_metrics[, "c_test"]

performance_table <- data.frame(
  Metric = c(
    "Apparent AUC",
    "Optimism-corrected AUC",
    "Brier score",
    "Optimism-corrected Brier score",
    "Apparent C-index",
    "Optimism-corrected C-index"
  ),
  Estimate = c(
    auc_app,
    auc_corrected,
    brier_app,
    brier_corrected,
    c_index_app,
    c_index_corrected
  ),
  CI_lower = c(
    auc_ci[1],
    quantile(boot_metrics[, "auc_test"], 0.025, na.rm = TRUE),
    quantile(brier_boot_values, 0.025, na.rm = TRUE),
    quantile(brier_boot_values, 0.025, na.rm = TRUE),
    auc_ci[1],
    quantile(c_boot_values, 0.025, na.rm = TRUE)
  ),
  CI_upper = c(
    auc_ci[3],
    quantile(boot_metrics[, "auc_test"], 0.975, na.rm = TRUE),
    quantile(brier_boot_values, 0.975, na.rm = TRUE),
    quantile(brier_boot_values, 0.975, na.rm = TRUE),
    auc_ci[3],
    quantile(c_boot_values, 0.975, na.rm = TRUE)
  )
)

write.csv(performance_table,
          file.path(out_dir, "model_performance_metrics_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")


## =========================
## 8. 校准指标
## =========================

cal_model <- glm(df_model$MCVT_bin ~ safe_logit(pred_prob), family = binomial())
calibration_intercept <- coef(cal_model)[1]
calibration_slope <- coef(cal_model)[2]

cal_in_large_model <- glm(df_model$MCVT_bin ~ offset(lp), family = binomial())
calibration_in_the_large <- coef(cal_in_large_model)[1]

bootstrap_calibration_slope <- function(data, indices) {
  d_boot <- data[indices, ]

  fit_b <- tryCatch(
    glm(MCVT_bin ~ Time + NLR + PLR + TG_linear + TG_rcs1 + TC + APTT,
        data = d_boot, family = binomial()),
    error = function(e) NULL
  )

  if (is.null(fit_b)) return(c(slope_boot = NA, slope_test = NA))

  p_boot <- as.numeric(predict(fit_b, newdata = d_boot, type = "response"))
  p_test <- as.numeric(predict(fit_b, newdata = data, type = "response"))

  slope_boot <- tryCatch(
    coef(glm(d_boot$MCVT_bin ~ safe_logit(p_boot), family = binomial()))[2],
    error = function(e) NA_real_
  )

  slope_test <- tryCatch(
    coef(glm(data$MCVT_bin ~ safe_logit(p_test), family = binomial()))[2],
    error = function(e) NA_real_
  )

  c(slope_boot = slope_boot, slope_test = slope_test)
}

set.seed(2025)
boot_cal <- t(apply(boot_indices, 2, function(idx) {
  bootstrap_calibration_slope(df_model, idx)
}))

slope_optimism <- mean(boot_cal[, "slope_boot"] - boot_cal[, "slope_test"], na.rm = TRUE)
bootstrap_corrected_calibration_slope <- calibration_slope - slope_optimism

corrected_slope_distribution <- calibration_slope -
  (boot_cal[, "slope_boot"] - boot_cal[, "slope_test"])

calibration_table <- data.frame(
  Metric = c(
    "Calibration intercept",
    "Calibration slope",
    "Calibration-in-the-large",
    "Bootstrap-corrected calibration slope",
    "Brier score"
  ),
  Estimate = c(
    calibration_intercept,
    calibration_slope,
    calibration_in_the_large,
    bootstrap_corrected_calibration_slope,
    brier_app
  ),
  CI_lower = c(
    NA,
    NA,
    NA,
    quantile(corrected_slope_distribution, 0.025, na.rm = TRUE),
    quantile(brier_boot_values, 0.025, na.rm = TRUE)
  ),
  CI_upper = c(
    NA,
    NA,
    NA,
    quantile(corrected_slope_distribution, 0.975, na.rm = TRUE),
    quantile(brier_boot_values, 0.975, na.rm = TRUE)
  )
)

write.csv(calibration_table,
          file.path(out_dir, "calibration_metrics_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")


## =========================
## 9. Hosmer-Lemeshow 检验
## =========================

hl <- ResourceSelection::hoslem.test(
  x = df_model$MCVT_bin,
  y = pred_prob,
  g = 10
)

hl_table <- data.frame(
  Hosmer_Lemeshow_chi_square = as.numeric(hl$statistic),
  df = as.numeric(hl$parameter),
  P_value = as.numeric(hl$p.value)
)

write.csv(hl_table,
          file.path(out_dir, "hosmer_lemeshow_test_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")

hl_groups <- data.frame(hl$observed, hl$expected)
write.csv(hl_groups,
          file.path(out_dir, "hosmer_lemeshow_groups_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")


## =========================
## 10. DCA 及不确定性分析
## =========================

net_benefit <- function(y, p, thresholds) {
  n <- length(y)
  prevalence <- mean(y)
  out <- lapply(thresholds, function(pt) {
    pred_pos <- p >= pt
    tp <- sum(pred_pos & y == 1)
    fp <- sum(pred_pos & y == 0)

    nb_model <- tp / n - fp / n * pt / (1 - pt)
    nb_all <- prevalence - (1 - prevalence) * pt / (1 - pt)
    nb_none <- 0

    data.frame(
      threshold = pt,
      model = nb_model,
      treat_all = nb_all,
      treat_none = nb_none
    )
  })
  dplyr::bind_rows(out)
}

thresholds <- seq(0.01, 0.99, by = 0.01)
dca_app <- net_benefit(df_model$MCVT_bin, pred_prob, thresholds)

set.seed(2025)
dca_boot_array <- array(NA_real_, dim = c(length(thresholds), B))

for (b in seq_len(B)) {
  idx <- sample(seq_len(nrow(df_model)), replace = TRUE)
  d_b <- df_model[idx, ]
  p_b <- pred_prob[idx]
  nb_b <- net_benefit(d_b$MCVT_bin, p_b, thresholds)
  dca_boot_array[, b] <- nb_b$model
}

dca_result <- dca_app
dca_result$model_CI_lower <- apply(dca_boot_array, 1, quantile, probs = 0.025, na.rm = TRUE)
dca_result$model_CI_upper <- apply(dca_boot_array, 1, quantile, probs = 0.975, na.rm = TRUE)

write.csv(dca_result,
          file.path(out_dir, "dca_net_benefit_with_ci_updated.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")


## =========================
## 11. 汇总输出到 Excel
## =========================

wb <- openxlsx::createWorkbook()

openxlsx::addWorksheet(wb, "Univariate")
openxlsx::writeData(wb, "Univariate", univariate_table)

openxlsx::addWorksheet(wb, "LASSO")
openxlsx::writeData(wb, "LASSO", lasso_coef_table)
openxlsx::writeData(wb, "LASSO", data.frame(lambda_min = lambda_min,
                                             lambda_1se = lambda_1se),
                    startRow = nrow(lasso_coef_table) + 3)

openxlsx::addWorksheet(wb, "RCS")
openxlsx::writeData(wb, "RCS", rcs_result)

openxlsx::addWorksheet(wb, "Model coefficients")
openxlsx::writeData(wb, "Model coefficients", coef_table)

openxlsx::addWorksheet(wb, "EPV shrinkage")
openxlsx::writeData(wb, "EPV shrinkage", epv_table)

openxlsx::addWorksheet(wb, "Performance")
openxlsx::writeData(wb, "Performance", performance_table)

openxlsx::addWorksheet(wb, "Calibration")
openxlsx::writeData(wb, "Calibration", calibration_table)

openxlsx::addWorksheet(wb, "Hosmer-Lemeshow")
openxlsx::writeData(wb, "Hosmer-Lemeshow", hl_table)
openxlsx::writeData(wb, "Hosmer-Lemeshow", hl_groups, startRow = 5)

openxlsx::addWorksheet(wb, "DCA")
openxlsx::writeData(wb, "DCA", dca_result)

openxlsx::saveWorkbook(
  wb,
  file.path(out_dir, "updated_statistical_results_R.xlsx"),
  overwrite = TRUE
)


## =========================
## 12. 控制台打印关键结果
## =========================

cat("\n===== LASSO =====\n")
cat("lambda.min =", lambda_min, "\n")
cat("lambda.1se =", lambda_1se, "\n")
cat("Selected variables at lambda.1se:\n")
print(selected_lasso)

cat("\n===== RCS =====\n")
print(rcs_result)

cat("\n===== Final model coefficients =====\n")
print(coef_table)

cat("\n===== EPV and shrinkage =====\n")
print(epv_table)

cat("\n===== Model performance =====\n")
print(performance_table)

cat("\n===== Calibration =====\n")
print(calibration_table)

cat("\n===== Hosmer-Lemeshow test =====\n")
print(hl_table)

cat("\nAll results have been saved to:\n")
cat(out_dir, "\n")