survMisc::comp计算p值错误：协方差矩阵问题及修复咨询

问题背景 #

我在分析一个包含407个样本的种子萌发时间数据集，变量包括：

• time：萌发天数（1-7天）
• status：事件状态（1=萌发，0=截尾）
• stress：处理组标识（1=施加胁迫，0=未施加）
  数据为右截尾且不满足比例风险（PH）假设。

用survival::survdiff()和survMisc::comp()做log-rank检验时，得到的p值不一致：前者给出p=0.03，后者给出p=0.052。排查后发现问题出在survMisc::COV()计算协方差矩阵错误，它算出的VAR(Oi-Ei)为77.83，而手动计算和survdiff的正确结果是62.2。

尝试把数据转成宽格式后，COV()能得到正确的协方差矩阵，但没法直接传入comp()使用。推测问题和数据集特性（比如时间点少、事件数分布）以及COV()的内部公式有关，想找直接修复comp()问题的方法。

可行解决方案 #

1. 手动修正comp()的协方差计算 #

comp()的核心统计量依赖COV()的结果，我们可以手动替换正确的方差值，重新计算卡方和p值：

library(survival)
library(survMisc)

# 先获取survdiff的正确结果
sdiff_res <- survdiff(Surv(time, status) ~ stress, data = your_data)
correct_var <- sdiff_res$var[1,1] # 两组比较下，取第一行第一列的方差值

# 调用原comp()获取基础结果
comp_res <- comp(Surv(time, status) ~ stress, data = your_data)

# 替换协方差矩阵并重新计算统计量
comp_res$cov <- matrix(correct_var, nrow = 1, ncol = 1) # 两组比较时是1x1矩阵
comp_res$chisq <- t(comp_res$oe) %*% solve(comp_res$cov) %*% comp_res$oe
comp_res$p <- pchisq(comp_res$chisq, df = nrow(comp_res$cov), lower.tail = FALSE)

# 查看修正后的结果
comp_res

2. 调整数据集结构适配COV() #

既然宽格式下COV()计算正确，可以把宽格式数据转成comp()能识别的长格式：

library(tidyr)

# 假设宽格式数据是每行一个样本，列包含time1到time7的状态（1=萌发，0=未萌发/截尾）
long_data <- pivot_longer(your_wide_data, 
                          cols = starts_with("time"), 
                          names_to = "time", 
                          values_to = "status")
# 提取时间数值
long_data$time <- as.integer(sub("time", "", long_data$time))
# 确保status定义正确：已萌发的样本在对应时间点标记为1，后续时间点设为0避免重复计数
long_data <- long_data %>% 
  group_by(id) %>% 
  mutate(status = ifelse(cumsum(status) > 0, 0, status)) %>% 
  ungroup()

# 现在调用comp()
comp(Surv(time, status) ~ stress, data = long_data)

3. 拆分功能，组合使用工具包 #

如果comp()主要是用来做生存曲线可视化或汇总，完全可以用survdiff的正确统计量搭配survMisc的绘图功能：

sf <- survfit(Surv(time, status) ~ stress, data = your_data)
# 用survMisc绘制生存曲线
plot(sf)
# 标注survdiff的统计结果
text(x = 3, y = 0.8, labels = paste("Log-rank p =", round(sdiff_res$pvalue, 3)))

问题根源推测 #

survMisc::COV()的内部方差计算公式，在处理时间点少（仅7天）、事件分布相对集中的右截尾数据时，求和逻辑没有正确匹配每个时间点的风险人数和事件数，导致方差被高估。转宽格式后数据结构的改变，让COV()的求和逻辑正确对应了数据的实际事件发生情况，因此得到了正确结果。

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