R语言生存结局Brier分数计算及pec包维度报错解决方案

• 核心原因1：pec()的object参数默认接收拟合完成的生存模型对象（如coxph、survreg的输出），你直接传入预计算的单个时间点生存概率矩阵，不符合参数要求。
• 核心原因2：pec()默认exact=FALSE，此时即使你指定了times=12，函数仍会自动纳入所有小于maxtime的唯一事件时间点计算误差，你只提供了1个时间点的预测值，维度不匹配触发报错。

分两种常见场景选择对应方案：

场景1：用pec完整建模流程计算Brier score #

如果不需要提前预计算12个月生存概率，直接把指标作为自变量拟合生存模型传入pec()即可：

# 拟合index1对应的Cox预测模型，必须设置x=TRUE、y=TRUE供pec调用
cox_index1 <- coxph(Surv(time, status) ~ index1, data = data, x = TRUE, y = TRUE)

# 计算12个月的Brier score，加exact=TRUE指定只计算目标时间点
PredError <- pec(object = cox_index1,
                 formula = Surv(time, status) ~ 1,
                 cens.model = "marginal",
                 data = data,
                 verbose = F,
                 times = 12,
                 exact = TRUE)

# 提取12个月的Brier score值
PredError$AppErr[[1]]

如果要同时评估多个指标，把多个模型放进列表传入object即可：

cox_index2 <- coxph(Surv(time, status) ~ index2, data = data, x = TRUE, y = TRUE)
PredError <- pec(object = list(index1 = cox_index1, index2 = cox_index2),
                 formula = Surv(time, status) ~ 1,
                 cens.model = "marginal",
                 data = data,
                 verbose = F,
                 times = 12,
                 exact = TRUE)

场景2：已预计算好12个月生存概率，不需要重新建模 #

自定义预测函数告诉pec直接读取你预计算的概率值即可：

# 自定义预测函数，直接返回对应列的预计算概率
predict_precalc <- function(object, newdata, times){
  # object为存储预计算概率的列名，返回矩阵维度为 样本数*时间点数
  matrix(newdata[[object]], ncol = length(times))
}

# 调用pec，传入列名、自定义预测函数和exact=TRUE
PredError <- pec(object = c("index1", "index2"),
                 formula = Surv(time, status) ~ 1,
                 cens.model = "marginal",
                 data = data,
                 verbose = F,
                 times = 12,
                 exact = TRUE,
                 predictfun = predict_precalc)

# 提取两个指标对应的12个月Brier score
PredError$AppErr$index1
PredError$AppErr$index2

如果你的数据无截尾、不需要做逆概率加权调整，也可以直接手动计算未调整的Brier score：

brier_index1 <- mean((data$status[data$time <=12] * (1 - data$index1[data$time <=12])^2) + ((1 - data$status[data$time <=12]) * (data$index1[data$time <=12])^2))

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