coxme包一致性与森林图实现或survival包多分组语法求助

我正在拟合迟滞进入生存模型（staggered entry survival models），个体存在差异化监测时间，模型需要按individualid分组，同时个体隶属于更大的组（比如同一鸟巢的多个个体），groupid需作为随机效应。

使用survival包建模时遇到多分组问题：无法同时使用frailty项与cluster项；添加两个frailty项时直接报错。转而用coxme包可以构建包含groupid和individualid的混合效应模型，但无法计算一致性（concordance）或生成森林图（ggforest）。

寻求以下两种解决方案之一：

• A）在survival包中实现多分组的方法
• B）在coxme包中计算一致性与生成森林图的方法

示例代码 #

library(survival)
library(survminer)
library(coxme)

# 创建数据集
individualid <- c(1,1,2,2,2,3,3,3,3,3,3,3,4,5,6,6,6)
groupid <- c(1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2)
event <- c(0,1,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0)
entry <- c(0,1,0,1,2,0,1,2,3,4,5,6,0,0,0,1,2)
exit <- entry + 1
precip <- c(1, 0, 1, 0, .2, 1, 0, 0.2, 0,0, 0, 0, 3, 3, 3, 2, 0)
temp <- c(50, 70, 50, 70, 62, 50, 70, 62, 44, 81, 92, 66, 32, 32, 32, 51, 42)

data <- data.frame(individualid, groupid, event, entry, exit, precip, temp)

# 单变量模型
model1 <- coxph(Surv(entry, exit, event) ~ temp,
                data = data)

# 带随机截距的模型（groupid）
model2 <- coxph(Surv(entry, exit, event) ~ temp +
                  frailty(groupid),
                data = data)

# 按individualid分组的模型
model3 <- coxph(Surv(entry, exit, event) ~ temp, 
                cluster = data$individualid,
                data = data)

# 尝试同时加入frailty和cluster的模型
model4 <- coxph(Surv(entry, exit, event) ~ temp + frailty(groupid), 
                cluster = data$individualid,
                data = data)
# 该模型在示例数据中运行，但真实数据报错，survival官方文档明确禁止frailty与cluster同时使用

# 尝试加入两个frailty项的模型
# model5 <- coxph(Surv(entry, exit, event) ~ temp + frailty(groupid) + frailty(individualid),
#                 data = data)
# 示例数据中运行会崩溃，真实数据报错：
# Error in coxpenal.fit(...) : Only one sparse penalty term allowed

# 模型1的一致性计算
concordance(model1)
# 正常运行

# 模型1的森林图
survminer::ggforest(model1, data = as.data.frame(data))
# 正常运行

# 使用coxme构建混合效应模型（groupid和individualid为随机效应）
model6 <- coxme(Surv(entry, exit, event) ~ temp + (1|groupid) + (1|individualid), 
                data = data)

model6$concordance
# 返回空值

survminer::ggforest(model6, data = as.data.frame(data))
# 报错

解决方案A：survival包实现多分组的替代方案 #

survival包的coxph确实限制了同时使用frailty和cluster，且仅支持单个frailty项。可通过以下两种方式处理两级分组：

1. 分层+Frailty组合 #

将individualid作为分层变量，强制模型在每个个体层内拟合，同时用frailty(groupid)处理组水平随机效应：

model_A <- coxph(Surv(entry, exit, event) ~ temp + frailty(groupid) + strata(individualid),
                 data = data)

注：该方法假设个体水平效应适合作为固定分层而非随机效应，若必须将个体设为随机效应，此方法不适用。

2. Frailty+稳健标准误组合 #

用frailty(groupid)处理组随机效应，手动计算基于individualid的稳健标准误校正个体聚类影响（需注意统计假设合理性）：

library(sandwich)
library(lmtest)

model_A2 <- coxph(Surv(entry, exit, event) ~ temp + frailty(groupid),
                  data = data)
# 基于individualid计算稳健标准误
coeftest(model_A2, vcov = vcovCL, cluster = ~individualid)

解决方案B：coxme模型的一致性计算与森林图生成 #

1. 计算一致性（Concordance） #

coxme无内置一致性方法，可手动计算或用第三方包：

方法1：手动计算

# 提取线性预测值
lp <- predict(model6, type = "linear.predictors")
# 构建生存对象
surv_obj <- with(data, Surv(entry, exit, event))
# 使用survival包的concordance函数
concordance(surv_obj ~ lp, data = data)

方法2：用pec包计算

library(pec)
cindex(model6, formula = Surv(entry, exit, event) ~ 1, data = data)

2. 生成森林图（ggforest） #

survminer::ggforest不支持coxme对象，可手动提取系数信息用ggplot2绘制：

library(ggplot2)

# 提取固定效应的系数、置信区间
coef_df <- as.data.frame(summary(model6)$coefficients)
coef_df$term <- rownames(coef_df)
coef_df$lower <- coef_df$`exp(coef)` - 1.96 * coef_df$`se(coef)`
coef_df$upper <- coef_df$`exp(coef)` + 1.96 * coef_df$`se(coef)`

# 绘制森林图
ggplot(coef_df, aes(x = term, y = `exp(coef)`)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
  coord_flip() +
  labs(title = "CoxME模型森林图",
       y = "风险比（95%置信区间）",
       x = "协变量") +
  theme_minimal()

若需展示随机效应方差，可额外提取summary(model6)$varcor的信息添加到图中。

内容的提问来源于stack exchange，提问作者blake.struthers