如何基于cph系数手动构建Cox PH模型并获取模型方程？

一、先明确Cox PH模型的核心方程 #

Cox比例风险模型的生存概率计算公式为：
$$S(t|X) = S_0(t)^{\exp(\beta_1 var1 + \beta_2 var2 + \beta_3 var3 + \beta_4 var4)}$$
其中：

• $S_0(t)$：基准生存概率（所有协变量取0时的生存概率）
• $\beta_1 \sim \beta_4$：模型的回归系数
• $var1 \sim var4$：新个体的协变量值

要实现无原始数据库的预测，只需从已训练好的cph模型中提取$\beta$系数和$S_0(t)$即可。

二、从已训练的cph模型提取关键参数 #

在原建模环境中执行以下代码，提取并保存所需参数：

library(rms)

# 假设cox1是已训练好的cph模型
# 1. 提取回归系数
beta_coef <- coef(cox1)
# 2. 提取指定时间点的基准生存概率（对应12/24/36个月）
base_surv <- survest(cox1, 
                     newdata = data.frame(var1=0, var2=0, var3=0, var4=0),
                     times = c(12,24,36))$surv
# 3. 保存参数到本地（供Shiny调用）
saveRDS(list(beta = beta_coef, 
             base_surv = base_surv, 
             target_times = c(12,24,36),
             var_names = names(beta_coef)), 
        "cox_model_params.rds")

三、Shiny应用中实现无依赖预测 #

这里提供两种实现方式，推荐第一种（手动计算，无需重构复杂模型对象）：

方式1：手动计算生存概率（最稳妥） #

library(shiny)

# 加载预存的模型参数
model_params <- readRDS("cox_model_params.rds")
beta <- model_params$beta
base_surv <- model_params$base_surv
target_times <- model_params$target_times
var_names <- model_params$var_names

ui <- fluidPage(
  titlePanel("Cox模型生存概率预测"),
  sidebarLayout(
    sidebarPanel(
      numericInput("var1", "var1值:", value = 90),
      numericInput("var2", "var2值:", value = 1),
      numericInput("var3", "var3值:", value = 8),
      numericInput("var4", "var4值:", value = 8)
    ),
    mainPanel(
      tableOutput("surv_result")
    )
  )
)

server <- function(input, output) {
  output$surv_result <- renderTable({
    # 构建协变量向量
    X <- c(input$var1, input$var2, input$var3, input$var4)
    names(X) <- var_names
    
    # 计算指数化线性预测值
    exp_linear <- exp(sum(X * beta))
    # 计算各时间点生存概率
    surv_probs <- base_surv ^ exp_linear
    
    # 整理输出表格
    data.frame(
      时间(月) = target_times,
      生存概率 = round(surv_probs, 4)
    )
  })
}

shinyApp(ui, server)

方式2：重构cph模型对象（复用predictSurvProb） #

如果想继续使用predictSurvProb函数，可以重构cph模型对象（需保留足够多的模型属性）：

# 原环境中先保存完整模型核心组件
cox_recon <- list(
  coefficients = coef(cox1),
  terms = cox1$terms,
  Design = cox1$Design,
  xlevels = cox1$xlevels,
  surv = cox1$surv
)
class(cox_recon) <- "cph"
saveRDS(cox_recon, "cox_reconstructed.rds")

# Shiny应用中调用
library(shiny)
library(pec)
library(rms)

cox_recon <- readRDS("cox_reconstructed.rds")

ui <- fluidPage(
  # 同方式1的UI
)

server <- function(input, output) {
  output$surv_result <- renderTable({
    new_patient <- data.frame(
      var1 = input$var1,
      var2 = input$var2,
      var3 = input$var3,
      var4 = input$var4
    )
    surv_probs <- predictSurvProb(cox_recon, new_patient, times = c(12,24,36))
    data.frame(
      时间(月) = c(12,24,36),
      生存概率 = round(t(surv_probs), 4)
    )
  })
}

shinyApp(ui, server)

注意事项 #

1. 如果原模型中存在因子型变量，需在保存参数时额外保存xlevels，并在Shiny中把输入转换为对应因子水平，否则预测结果会出错。
2. 若基准生存函数未覆盖目标时间点，survest会自动插值计算对应概率，无需额外处理。

内容的提问来源于stack exchange，提问作者B_slash_