混合效应Cox回归及时间交互项的校正生存曲线绘制问询

Q1：coxme模型固定效应校正生存曲线绘制 #

对于带随机效应的coxme模型，survminer等工具确实不直接支持绘制固定效应校正曲线，但可以手动提取固定效应系数、基线风险，进而计算生存概率并绘图，具体实现如下：

代码示例 #

library(survival)
library(coxme)
library(ggplot2)

# 拟合coxme模型
kidney <- data.frame(kidney)
coxme_mod <- coxme(Surv(time, status) ~ age + sex + (1|disease), data = kidney)

# 提取固定效应系数
fixed_coef <- fixef(coxme_mod)

# 构造参考数据集：age取均值，sex覆盖所有水平（假设sex为0/1二分类）
ref_data <- expand.grid(
  age = mean(kidney$age, na.rm = TRUE),
  sex = c(0, 1),
  disease = unique(kidney$disease)[1] # 随机效应后续会排除，取值不影响结果
)

# 生成设计矩阵并计算线性预测值（排除随机效应）
X <- model.matrix(~ age + sex, data = ref_data)
lp <- as.vector(X %*% fixed_coef)

# 获取群体平均基线累积风险（coxme已自动积分掉随机效应）
base_haz <- basehaz(coxme_mod, centered = FALSE)
colnames(base_haz) <- c("time", "hazard")

# 计算各参考组的生存概率
surv_df <- lapply(seq_along(lp), function(i) {
  data.frame(
    time = base_haz$time,
    survival = exp(-exp(lp[i]) * base_haz$hazard),
    sex = factor(ref_data$sex[i], labels = c("Female", "Male"))
  )
})
surv_df <- do.call(rbind, surv_df)

# 绘制校正生存曲线
ggplot(surv_df, aes(x = time, y = survival, color = sex)) +
  geom_line(linewidth = 1) +
  labs(x = "时间", y = "生存概率", color = "性别") +
  theme_bw()

关键说明 #

• fixef(coxme_mod)专门提取模型的固定效应系数，完全排除随机效应的干扰；
• basehaz(coxme_mod)返回的是群体平均水平的基线累积风险，已经自动处理了随机效应的积分，适合用于绘制固定效应校正曲线；
• 生存概率计算公式遵循Cox模型定义：S(t) = exp(-exp(线性预测值) × 基线累积风险(t))。

Q2：含时间交互项+frailty的coxph模型校正曲线绘制 #

带tt()时间交互项的模型，线性预测值随时间动态变化，同样可以手动计算生存概率并绘图，具体实现如下：

代码示例 #

library(survival)
library(ggplot2)

# 拟合带时间交互和frailty的coxph模型
kidney <- data.frame(kidney)
coxph_tt_mod <- coxph(
  Surv(time, status) ~ age + sex + tt(sex) + frailty(disease), 
  data = kidney,
  tt=function(x,t,...) x*t
)

# 提取模型系数
coefs <- coef(coxph_tt_mod)

# 构造参考数据集：age取均值，sex覆盖所有水平
ref_data <- expand.grid(
  age = mean(kidney$age, na.rm = TRUE),
  sex = c(0, 1)
)

# 获取基线累积风险（frailty设为0时的基线风险）
base_haz <- basehaz(coxph_tt_mod, centered = FALSE)
colnames(base_haz) <- c("time", "hazard")

# 计算各参考组的生存概率（线性预测值随时间变化）
surv_df <- lapply(1:nrow(ref_data), function(i) {
  age_val <- ref_data$age[i]
  sex_val <- ref_data$sex[i]
  # 每个时间点的线性预测值：包含sex与时间的交互项
  lp_t <- coefs["age"]*age_val + coefs["sex"]*sex_val + coefs["tt(sex)"]*sex_val*base_haz$time
  # 计算生存概率
  survival <- exp(-exp(lp_t) * base_haz$hazard)
  data.frame(
    time = base_haz$time,
    survival = survival,
    sex = factor(sex_val, labels = c("Female", "Male"))
  )
})
surv_df <- do.call(rbind, surv_df)

# 绘制生存曲线
ggplot(surv_df, aes(x = time, y = survival, color = sex)) +
  geom_line(linewidth = 1) +
  labs(x = "时间", y = "生存概率", color = "性别") +
  theme_bw()

关键说明 #

• tt(sex)定义为sex*t，因此线性预测值随时间t动态变化，需要针对每个时间点重新计算；
• frailty(disease)作为随机效应，计算时默认排除（即设为0，对应exp(0)=1，不改变风险比），得到的是固定效应校正后的生存曲线；
• 若需要考虑frailty的影响，可以在计算时加入frailty的预测值，但校正曲线通常仅关注固定效应的组间对比。

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