#' Cox Proportional Hazards Regression
#'
#' @export
coxp <- function(dataset,
                 time,
                 status,
                 evar,
                 int = "",
                 check = "",
                 form,
                 data_filter = "",
                 arr = "",
                 rows = NULL,
                 envir = parent.frame()) {
  
  if (!requireNamespace("survival", quietly = TRUE))
    stop("survival package is required but not installed.")
  attachNamespace("survival") 
  on.exit(detach("package:survival"), add = TRUE)
  
  ## ---- 公式入口 ----------------------------------------------------------
  if (!missing(form)) {
    form  <- as.formula(format(form))
    vars  <- all.vars(form)
    time  <- vars[1]
    status<- vars[2]
    evar  <- vars[-(1:2)]
  }
  
  ## ---- 基础检查 ----------------------------------------------------------
  if (time %in% evar || status %in% evar) {
    return("Time/status variable contained in explanatory variables." %>%
             add_class("coxp"))
  }
  
  vars <- unique(c(time, status, evar))
  df_name <- if (is_string(dataset)) dataset else deparse(substitute(dataset))
  dataset <- get_data(dataset, vars, filt = data_filter, arr = arr, rows = rows, envir = envir)
  
  ## 状态变量检查与转换
  surv_status <- dataset[[status]]
  if (!is.numeric(surv_status)) {
    ## 允许 0/1、FALSE/TRUE、factor(未事件/事件) 等常见编码
    if (is.factor(surv_status) || is.character(surv_status)) {
      lv <- unique(surv_status)
      if (length(lv) != 2) {
        return("Status variable must be binary (0/1 or two levels)." %>% add_class("coxp"))
      }
      dataset[[status]] <- as.numeric(factor(surv_status, levels = lv)) - 1L
    } else {
      return("Status variable must be numeric 0/1 or binary factor." %>% add_class("coxp"))
    }
  } else {
    if (!all(unique(surv_status) %in% c(0, 1))) {
      return("Status variable must contain only 0 and 1." %>% add_class("coxp"))
    }
  }
  
  if (missing(form)) {
    rhs <- if (length(evar) == 0) "1" else paste(evar, collapse = " + ")
    if (!is.empty(int)) rhs <- paste(rhs, paste(int, collapse = " + "), sep = " + ")
    form <- as.formula(paste("Surv(", time, ", ", status, ") ~ ", rhs))
  }
  
  if ("robust" %in% check) {
    model <- survival::coxph(form, data = dataset, robust = TRUE, x = TRUE, y = TRUE)
  } else {
    model <- survival::coxph(form, data = dataset, x = TRUE, y = TRUE)
  }
  
  ## 失败模型保护
  if (inherits(model, "try-error")) {
    return("Model estimation failed. Check data separation or collinearity." %>% add_class("coxp"))
  }
  
  ## 基础摘要信息
  coef_df <- broom::tidy(model, conf.int = TRUE)  # 系数、HR、CI、p
  n  <- nrow(dataset)          # 样本量
  n_event <- sum(dataset[[status]])  # 事件数
  conc <- tryCatch(
    survival::concordancefit(
      y = Surv(dataset[[time]], dataset[[status]]),
      x = predict(model, type = "lp"),
      data = dataset
    )$concordance,
    error = function(e) NA
  )
  ## 打包返回
  out <- as.list(environment())
  out$model   <- model
  out$df_name <- df_name
  out$type    <- "survival"
  out$check   <- check
  ## 附加对象
  out$coef_df <- coef_df
  out$n       <- n
  out$n_event <- n_event
  out$concordance <- conc
  out$model$model <- dataset[, all.vars(form), drop = FALSE]
  out$time_var  <- time
  out$status_var <- status
  add_class(out, c("coxp", "model"))
}

#' @export
summary.coxp <- function(object, dec = 3, ...) {
  if (is.character(object)) {
    cat(object, "\n")
    return(invisible(object))
  }
  
  if (!inherits(object$model, "coxph")) {
    cat("** Invalid Cox model object. **\n")
    return(invisible(object))
  }
  
  ## 基础模型信息
  cat("Cox Proportional Hazards Regression\n")
  cat("Data:", object$df_name, "\n")
  if (!is.empty(object$data_filter)) {
    cat("Filter:", gsub("\\n", "", object$data_filter), "\n")
  }
  if (!is.empty(object$arr)) {
    cat("Arrange:", gsub("\\n", "", object$arr), "\n")
  }
  if (!is.empty(object$rows)) {
    cat("Slice:", gsub("\\n", "", object$rows), "\n")
  }
  cat("Time variable   :", object$time, "\n")
  cat("Status variable :", object$status, "\n")
  cat("Explanatory vars:", paste(object$evar, collapse = ", "), "\n")
  cat("N =", object$n, ", Events =", object$n_event, "\n")
  cat("Concordance =", sprintf("%.3f", object$concordance), "\n\n")
  
  ## 系数表
  coef_df <- object$coef_df
  coef_df$sig_star <- sig_stars(coef_df$p.value) %>% format(justify = "left")
  coef_df$label <- rownames(coef_df)
  
  ## 格式化输出
  coeff <- coef_df %>%
    mutate(
      HR = sprintf("%.3f", exp(estimate)),
      `HR.low` = sprintf("%.3f", exp(conf.low)),
      `HR.high` = sprintf("%.3f", exp(conf.high)),
      coef = sprintf("%.3f", estimate),
      se = sprintf("%.3f", std.error),
      z = sprintf("%.3f", statistic),
      p = ifelse(p.value < .001, "< .001", sprintf("%.3f", p.value))
    ) %>%
    select(label, coef, se, z, p, sig_star, HR, HR.low, HR.high)
  
  colnames(coeff) <- c(" ", "Coef", "SE", "z", "p", " ", "HR", "HR.lower", "HR.upper")
  print.data.frame(coeff, row.names = FALSE)
  
  cat("\nSignif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n")
  
  ## 模型检验
  sm <- summary(object$model)
  cat("Likelihood ratio test =", sprintf("%.2f", sm$logtest[1]), "on", sm$logtest[2], "df, p =",
      ifelse(sm$logtest[3] < .0001, "< .0001", sprintf("%.4f", sm$logtest[3])), "\n")
  cat("Wald test             =", sprintf("%.2f", sm$waldtest[1]), "on", sm$waldtest[2], "df, p =",
      ifelse(sm$waldtest[3] < .0001, "< .0001", sprintf("%.4f", sm$waldtest[3])), "\n")
  cat("Score (logrank) test  =", sprintf("%.2f", sm$sctest[1]), "on", sm$sctest[2], "df, p =",
      ifelse(sm$sctest[3] < .0001, "< .0001", sprintf("%.4f", sm$sctest[3])), "\n")
  
  invisible(object)
}


#' @export
predict.coxp <- function(object, pred_data = NULL, pred_cmd = "",
                         conf_lev = 0.95, dec = 3, envir = parent.frame(), ...) {
  if (is.character(object)) return(object)
  
  # 1. 构造预测数据
  if (is.null(pred_data)) {
    newdata <- envir$.model_frame %||% object$model$model
  } else {
    # 获取预测数据集
    newdata <- get_data(pred_data, vars = NULL, envir = envir)
    
    # 变量存在性校验
    model_evar <- object$evar  
    pred_cols  <- colnames(newdata)
    missing_vars <- setdiff(model_evar, pred_cols)
    
    if (length(missing_vars) > 0) {
      return(paste0(
        "All variables in the model must also be in the prediction data\n",
        "Variables in the model: ", paste(model_evar, collapse = ", "), "\n",
        "Variables not available in prediction data: ", paste(missing_vars, collapse = ", ")
      ) %>% add_class("coxp.predict"))
    }
    newdata <- newdata[, model_evar, drop = FALSE]
  }
  
  # 2. 应用 pred_cmd
  if (!is.empty(pred_cmd)) {
    newdata <- modify_data(newdata, pred_cmd, envir = envir)
  }
  
  # 3. 过滤全NA的列
  newdata <- newdata[, colSums(is.na(newdata)) != nrow(newdata), drop = FALSE]
  if (ncol(newdata) == 0 && length(object$evar) > 0) {
    return(paste("预测数据中无有效解释变量列（需包含：", paste(object$evar, collapse = ", "), "）") %>% add_class("coxp.predict"))
  }
  
  # 4. 核心预测计算
  pred_cox <- predict(
    object$model, 
    newdata = newdata, 
    type = "lp", 
    se.fit = TRUE
  )
  z_val <- qnorm((1 + conf_lev) / 2)
  lp_lower <- pred_cox$fit - z_val * pred_cox$se.fit
  lp_upper <- pred_cox$fit + z_val * pred_cox$se.fit
  hr <- exp(pred_cox$fit)
  hr_lower <- exp(lp_lower)
  hr_upper <- exp(lp_upper)
  
  # 5. 构建结果数据框
  pred_result <- data.frame(
    lp = round(pred_cox$fit, dec),
    HR = round(hr, dec),
    lp_lower = round(lp_lower, dec),
    lp_upper = round(lp_upper, dec),
    HR_lower = round(hr_lower, dec),
    HR_upper = round(hr_upper, dec),
    check.names = FALSE,
    stringsAsFactors = FALSE
  )
  pred_full <- cbind(newdata, pred_result)
  
  # 6. 添加元信息
  pred_full <- pred_full %>%
    radiant.data::set_attr("radiant_df_name", object$df_name) %>%
    radiant.data::set_attr("radiant_time", object$time) %>%
    radiant.data::set_attr("radiant_status", object$status) %>%
    radiant.data::set_attr("radiant_evar_actual", colnames(newdata)) %>%
    radiant.data::set_attr("radiant_conf_lev", conf_lev) %>%
    radiant.data::set_attr("radiant_dec", dec) %>%
    add_class(c("coxp.predict", "model.predict"))
  
  return(pred_full)
}


#' @export
print.coxp.predict <- function(x, ..., n = 10) {
  if (is.character(x)) {
    cat(x, "\n")
    return(invisible(x))
  }
  
  # 转为数据框
  x_df <- as.data.frame(x, stringsAsFactors = FALSE)
  
  df_name <- attr(x_df, "radiant_df_name") %||% "Unknown"
  time_var <- attr(x_df, "radiant_time") %||% "Unknown"
  status_var <- attr(x_df, "radiant_status") %||% "Unknown"
  conf_lev <- attr(x_df, "radiant_conf_lev") %||% 0.95
  dec <- attr(x_df, "radiant_dec") %||% 3
  ci_perc <- paste0(c(round((1 - conf_lev) / 2 * 100, 1), 
                      round((1 + conf_lev) / 2 * 100, 1)), "%")
  
  total_cols <- ncol(x_df)
  result_count <- 6 
  
  if (total_cols < result_count) {
    cat("Error: Not enough columns for prediction results (need at least 6 result columns).\n")
    return(invisible(x))
  }
  
  evar_count_actual <- total_cols - result_count
  evar_cols_actual <- colnames(x_df)[1:evar_count_actual]
  
  new_result_names <- c(
    "lp",
    "HR",
    ci_perc[1],                    
    ci_perc[2],                    
    paste0("HR_", ci_perc[1]),     
    paste0("HR_", ci_perc[2])      
  )
  
  # 拼接完整列名向量
  new_colnames <- c(evar_cols_actual, new_result_names)
  
  # 最终校验
  if (length(new_colnames) != total_cols) {
    cat("Error: Column name length mismatch.\n")
    cat("Total columns:", total_cols, "\n")
    cat("Constructed names:", length(new_colnames), "\n")
    cat("evar_cols_actual:", paste(evar_cols_actual, collapse = ", "), "\n")
    cat("new_result_names:", paste(new_result_names, collapse = ", "), "\n")
    return(invisible(x))
  }
  
  colnames(x_df) <- new_colnames
  
  cat("Cox Proportional Hazards Regression\n")
  cat("Data                 :", df_name, "\n")
  cat("Time variable        :", time_var, "\n")
  cat("Status variable      :", status_var, "\n")
  cat("Explanatory variables:", if (length(evar_cols_actual) > 0) paste(evar_cols_actual, collapse = ", ") else "None", "\n")
  cat("Confidence level     :", paste0(conf_lev * 100, "%"), "\n")
  cat("Total columns        :", total_cols, "(Explanatory:", evar_count_actual, ", Result:", result_count, ")\n")
  
  total_rows <- nrow(x_df)
  if (n == -1 || total_rows <= n) {
    cat("Rows shown           :", total_rows, "of", total_rows, "\n")
    out_df <- x_df
  } else {
    cat("Rows shown           :", n, "of", total_rows, "\n")
    out_df <- head(x_df, n)
  }
  cat("\n")
  
  # 格式化数值列
  numeric_cols <- (evar_count_actual + 1):total_cols
  if (length(numeric_cols) > 0) {
    out_df[, numeric_cols] <- lapply(out_df[, numeric_cols, drop = FALSE], function(col) {
      sprintf(paste0("%.", dec, "f"), as.numeric(col))
    })
  }
  
  print(out_df, row.names = FALSE)
  invisible(x)
}

#' @export
store.coxp.predict <- function(dataset, object, name = "hr", ...) {
  if (is.empty(name)) name <- "hr"
  name <- unlist(strsplit(name, "(\\s*,\\s*|\\s*;\\s*|\\s+)")) %>% 
    gsub("\\s", "", .) %>% 
    .[1] 
  
  pred_col <- "HR"
  if (!pred_col %in% colnames(object)) {
    stop("Prediction column 'HR' not found in prediction object.")
  }
  
  pred_df <- object[, pred_col, drop = FALSE]
  colnames(pred_df) <- name
  
  evar_actual <- attr(object, "radiant_evar_actual") %||% 
    attr(object, "radiant_evar") %||% character(0)
  
  indr <- indexr(dataset, vars = evar_actual, filt = "", rows = NULL, 
                 cmd = attr(object, "radiant_pred_cmd"))
  
  out_df <- as.data.frame(matrix(NA, nrow = nrow(dataset), ncol = 1), stringsAsFactors = FALSE)
  out_df[indr$ind, 1] <- pred_df[[1]]
  colnames(out_df) <- name
  
  dataset[, name] <- out_df
  dataset
}

#' @export
plot.coxp <- function(x, plots = "none", incl = NULL, incl_int = NULL,
                      conf_lev = 0.95, intercept = FALSE, lines = "",
                      nrobs = -1, shiny = FALSE, custom = FALSE, ...) {
  # 输入验证：检查模型是否有效
  if (is.character(x)) return(x)
  if (is.null(x$model)) {
    return("** 无效的模型输入：未找到模型数据 **")
  }
  
  # 检查是否有选择的图表类型
  if (length(plots) == 0 || all(plots == "none")) {
    if (shiny) {
      return("** 无图表可显示。请选择图表类型或重新运行计算 **")
    } else {
      return(invisible())
    }
  }
  
  plot_list <- list()
  nrCol <- 2  # 默认2列布局
  
  # 1. 系数图 (coef)
  if ("coef" %in% plots) {
    nrCol <- 1
    tryCatch({
      coef_df <- broom::tidy(x$model, conf.int = TRUE, conf.level = conf_lev)
      coef_df$hr <- exp(coef_df$estimate)
      coef_df$hr_low <- exp(coef_df$conf.low)
      coef_df$hr_high <- exp(coef_df$conf.high)
      
      # 处理包含的变量
      if (length(incl) > 0) {
        incl_regex <- paste0("^(", paste(incl, collapse = "|"), ")")
        coef_df <- coef_df[grepl(incl_regex, coef_df$term), ]
      }
      
      if (nrow(coef_df) == 0) {
        plot_list[["coef"]] <- "** 无系数可绘制（可能已排除所有变量）**"
      } else {
        p <- ggplot2::ggplot(coef_df, 
                             ggplot2::aes(x = .data$term, y = .data$hr, 
                                          ymin = .data$hr_low, ymax = .data$hr_high)) +
          ggplot2::geom_pointrange() +
          ggplot2::geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
          ggplot2::scale_x_discrete(limits = rev(coef_df$term)) +
          ggplot2::coord_flip() +
          ggplot2::labs(x = "", y = "风险比 (HR)", 
                        title = "Cox模型系数（风险比）") +
          ggplot2::theme(axis.text.y = ggplot2::element_text(hjust = 0))
        plot_list[["coef"]] <- p
      }
    }, error = function(e) {
      plot_list[["coef"]] <- paste0("** 系数图生成失败：", e$message, " **")
    })
  }
  
  # 2. 分布直方图 (dist)
  if ("dist" %in% plots) {
    tryCatch({
      dist_data <- x$model$model
      if (is.null(x$evar) || length(x$evar) == 0) {
        plot_list[["dist"]] <- "** 无解释变量可绘制分布 **"
      } else {
        for (v in x$evar) {
          if (v %in% colnames(dist_data)) {
            p <- visualize(dist_data, xvar = v, bins = 10, custom = TRUE) + 
              ggplot2::labs(title = paste(v, "的分布"))
            plot_list[[paste0("dist_", v)]] <- p
          } else {
            plot_list[[paste0("dist_", v)]] <- paste0("** 变量 ", v, " 不在数据中 **")
          }
        }
      }
    }, error = function(e) {
      plot_list[["dist_error"]] <- paste0("** 分布图生成失败：", e$message, " **")
    })
  }
  
  # 3. 特征重要性图 (vip)
  if ("vip" %in% plots) {
    nrCol <- 1
    tryCatch({
      if (length(x$evar) < 2) {
        plot_list[["vip"]] <- "** 至少需要2个变量才能生成特征重要性图 **"
      } else {
        coef_df <- broom::tidy(x$model)
        if (!"statistic" %in% colnames(coef_df)) {
          plot_list[["vip"]] <- "** 无法获取统计量，无法生成特征重要性图 **"
        } else {
          coef_df$Importance <- abs(coef_df$statistic)  # 使用z统计量
          coef_df <- coef_df[order(coef_df$Importance, decreasing = TRUE), ]
          
          p <- visualize(coef_df, xvar = "term", yvar = "Importance", 
                         type = "bar", custom = TRUE) +
            ggplot2::coord_flip() + 
            ggplot2::labs(title = "特征重要性", x = "", y = "重要性") +
            ggplot2::theme(axis.text.y = ggplot2::element_text(hjust = 0))
          plot_list[["vip"]] <- p
        }
      }
    }, error = function(e) {
      plot_list[["vip"]] <- paste0("** 特征重要性图生成失败：", e$message, " **")
    })
  }
  
  # 4. 生存曲线图 (survival)
  if ("survival" %in% plots) {
    nrCol <- 1
    tryCatch({
      # 检查时间和事件变量是否存在
      if (is.null(x$time_var) || is.null(x$status_var)) {
        plot_list[["survival"]] <- "** 模型缺少时间变量或事件变量，无法生成生存曲线 **"
      } else if (!x$time_var %in% colnames(x$model$model)) {
        plot_list[["survival"]] <- paste0("** 时间变量 '", x$time_var, "' 不在模型数据中 **")
      } else if (!x$status_var %in% colnames(x$model$model)) {
        plot_list[["survival"]] <- paste0("** 事件变量 '", x$status_var, "' 不在模型数据中 **")
      } else {
        # 检查是否有事件发生
        event_count <- sum(x$model$model[[x$status_var]] == 1, na.rm = TRUE)
        if (event_count == 0) {
          plot_list[["survival"]] <- "** 数据中无事件发生（status全为0），无法生成生存曲线 **"
        } else {
          # 计算生存曲线
          surv_fit <- survival::survfit(x$model, conf.int = TRUE, conf.lev = conf_lev)
          if (is.null(surv_fit) || length(surv_fit$time) == 0) {
            plot_list[["survival"]] <- "** 无法计算生存曲线，请检查模型和数据 **"
          } else {
            # 构建生存数据框
            surv_df <- data.frame(
              time = surv_fit$time,
              surv = surv_fit$surv,
              lower = ifelse(exists("lower", surv_fit), surv_fit$lower, NA),
              upper = ifelse(exists("upper", surv_fit), surv_fit$upper, NA)
            )
            
            p <- ggplot2::ggplot(surv_df, ggplot2::aes(x = .data$time, y = .data$surv)) +
              ggplot2::geom_line(color = "blue", linewidth = 1) +
              ggplot2::geom_ribbon(ggplot2::aes(ymin = .data$lower, ymax = .data$upper),
                                   alpha = 0.2, fill = "blue", na.rm = TRUE) +
              ggplot2::labs(title = "基线生存函数",
                            x = "时间", y = "生存概率") +
              ggplot2::ylim(0, 1)
            plot_list[["survival"]] <- p
          }
        }
      }
    }, error = function(e) {
      plot_list[["survival"]] <- paste0("** 生存曲线生成失败：", e$message, " **")
    })
  }
  
  # 5. 累积风险图 (cumhaz)
  if ("cumhaz" %in% plots) {
    nrCol <- 1
    tryCatch({
      # 检查时间和事件变量是否存在
      if (is.null(x$time_var) || is.null(x$status_var)) {
        plot_list[["cumhaz"]] <- "** 模型缺少时间变量或事件变量，无法生成累积风险曲线 **"
      } else if (!x$time_var %in% colnames(x$model$model)) {
        plot_list[["cumhaz"]] <- paste0("** 时间变量 '", x$time_var, "' 不在模型数据中 **")
      } else if (!x$status_var %in% colnames(x$model$model)) {
        plot_list[["cumhaz"]] <- paste0("** 事件变量 '", x$status_var, "' 不在模型数据中 **")
      } else {
        # 检查是否有事件发生
        event_count <- sum(x$model$model[[x$status_var]] == 1, na.rm = TRUE)
        if (event_count == 0) {
          plot_list[["cumhaz"]] <- "** 数据中无事件发生（status全为0），无法生成累积风险曲线 **"
        } else {
          # 计算生存曲线用于转换为累积风险
          surv_fit <- survival::survfit(x$model)
          if (is.null(surv_fit) || length(surv_fit$time) == 0) {
            plot_list[["cumhaz"]] <- "** 无法计算累积风险，请检查模型和数据 **"
          } else {
            # 构建累积风险数据框，处理surv=0的情况
            cumhaz_df <- data.frame(
              time = surv_fit$time,
              cumhaz = ifelse(surv_fit$surv > 0 & !is.na(surv_fit$surv), 
                              -log(surv_fit$surv), NA)
            )
            
            # 检查是否有有效数据
            valid_rows <- sum(!is.na(cumhaz_df$cumhaz))
            if (valid_rows < 2) {
              plot_list[["cumhaz"]] <- "** 有效数据不足，无法生成累积风险曲线 **"
            } else {
              p <- ggplot2::ggplot(cumhaz_df, ggplot2::aes(x = .data$time, y = .data$cumhaz)) +
                ggplot2::geom_line(color = "red", linewidth = 1, na.rm = TRUE) +
                ggplot2::labs(title = "基线累积风险函数",
                              x = "时间", y = "累积风险")
              plot_list[["cumhaz"]] <- p
            }
          }
        }
      }
    }, error = function(e) {
      plot_list[["cumhaz"]] <- paste0("** 累积风险曲线生成失败：", e$message, " **")
    })
  }
  
  # 6. Schoenfeld残差图 (schoenfeld)
  if ("schoenfeld" %in% plots) {
    nrCol <- 1
    tryCatch({
      # 检查事件变量是否存在
      if (is.null(x$status_var)) {
        plot_list[["schoenfeld"]] <- "** 模型缺少事件变量，无法计算Schoenfeld残差 **"
      } else if (!x$status_var %in% colnames(x$model$model)) {
        plot_list[["schoenfeld"]] <- paste0("** 事件变量 '", x$status_var, "' 不在模型数据中 **")
      } else {
        # 检查事件数是否足够（至少2个）
        event_count <- sum(x$model$model[[x$status_var]] == 1, na.rm = TRUE)
        if (event_count < 2) {
          plot_list[["schoenfeld"]] <- "** 事件数不足（至少需要2个事件），无法检验PH假设 **"
        } else {
          # 计算Schoenfeld残差
          schoenfeld <- survival::cox.zph(x$model)
          if (is.null(schoenfeld) || is.null(schoenfeld$time) || is.null(schoenfeld$y)) {
            plot_list[["schoenfeld"]] <- "** 无法计算Schoenfeld残差（可能PH假设严重违反）**"
          } else {
            # 处理单变量和多变量情况
            residual_vec <- if (is.matrix(schoenfeld$y)) {
              schoenfeld$y[, 1]  # 多变量时取第一列
            } else {
              schoenfeld$y       # 单变量时直接使用向量
            }
            
            # 构建残差数据框
            resid_df <- data.frame(
              time = schoenfeld$time,
              residual = residual_vec
            )
            
            # 检查有效数据
            if (nrow(resid_df) < 2 || all(is.na(resid_df$residual))) {
              plot_list[["schoenfeld"]] <- "** 有效残差不足，无法生成Schoenfeld图 **"
            } else {
              p <- ggplot2::ggplot(resid_df, 
                                   ggplot2::aes(x = .data$time, y = .data$residual)) +
                ggplot2::geom_point(alpha = 0.6, na.rm = TRUE) +
                ggplot2::geom_smooth(method = "loess", color = "red", se = TRUE, na.rm = TRUE) +
                ggplot2::geom_hline(yintercept = 0, linetype = "dashed") +
                ggplot2::labs(title = "Schoenfeld残差 - 检验PH假设",
                              x = "时间", y = "标准化Schoenfeld残差")
              plot_list[["schoenfeld"]] <- p
            }
          }
        }
      }
    }, error = function(e) {
      plot_list[["schoenfeld"]] <- paste0("** Schoenfeld残差图生成失败：", e$message, " **")
    })
  }
  
  # 7. 鞅残差图 (martingale)
  if ("martingale" %in% plots) {
    nrCol <- 1
    tryCatch({
      # 计算鞅残差和线性预测值
      martingale_resid <- residuals(x$model, type = "martingale")
      linear_pred      <- predict(x$model, type = "lp")
      if (is.null(martingale_resid) || is.null(linear_pred)) {
        plot_list[["martingale"]] <- "** 模型无法计算残差或线性预测值（可能拟合失败）**"
      } else if (length(martingale_resid) != length(linear_pred)) {
        plot_list[["martingale"]] <- "** 残差与线性预测值长度不匹配，无法绘图 **"
      } else {
        resid_df <- data.frame(linear_pred = linear_pred,
                               martingale  = martingale_resid) %>%
          stats::na.omit()
        if (nrow(resid_df) < 10) {
          plot_list[["martingale"]] <- paste0("** 有效数据点不足（仅", nrow(resid_df), "个），无法生成有意义的图 **")
        } else {
          p <- ggplot2::ggplot(resid_df,
                               ggplot2::aes(x = .data$linear_pred, y = .data$martingale)) +
            ggplot2::geom_point(alpha = 0.6) +
            ggplot2::labs(title = "鞅残差 vs 线性预测值",
                          x = "线性预测值", y = "鞅残差")
          
          if ("line" %in% lines)
            p <- p + ggplot2::geom_smooth(method = "lm", se = FALSE,
                                          color = "red", linewidth = 1.2)
          
          if ("loess" %in% lines)
            p <- p + ggplot2::geom_smooth(method = "loess", color = "blue", se = TRUE)
          
          if ("jitter" %in% lines)
            p <- p + ggplot2::geom_jitter(width = 0.2, height = 0.2,
                                          alpha = 0.3, color = "darkgreen")
          
          plot_list[["martingale"]] <- p
        }
      }
    }, error = function(e) {
      plot_list[["martingale"]] <- paste0("** 鞅残差图生成失败：", e$message, " **")
    })
  }
  
  # 检查是否有可显示的图
  if (length(plot_list) == 0) {
    return("** 无法生成任何图表。请检查模型设置和数据有效性 **")
  }
  
  # 输出图表
  if (custom) {
    if (length(plot_list) == 1) plot_list[[1]] else plot_list
  } else {
    if (length(plot_list) == 1) {
      plot_list[[1]]
    } else {
      patchwork::wrap_plots(plot_list, ncol = nrCol)
    }
  }
}