rms包cph模型仅含交互项报错及可视化方法求助

问题原因 #

rms::cph在解析包含限制性立方样条（rcs）的单独交互项（:形式）时，会因内部变量名匹配逻辑问题触发报错。而使用*时包含了主效应，变量名解析路径正常，因此模型可运行。

一、让仅含交互项的模型正常运行 #

方法1：手动构造交互项（适配cph） #

提前将样条转换与交互项计算为新变量，规避公式解析问题：

# 计算alk.phos的限制性立方样条转换值
pbc2$rcs_alkphos <- rcs(pbc2$alk.phos, 4)
# 构造sex与样条项的交互项
pbc2$sex_x_rcs_alkphos <- pbc2$sex * pbc2$rcs_alkphos
# 构建模型
mod_adj <- cph(Surv(tstart, tstop, death==2) ~ sex + sex_x_rcs_alkphos + albumin + bili, 
               id = factor(id), x=TRUE, data = pbc2)

方法2：改用survival::coxph模型 #

coxph对单独交互项的兼容性更强，直接使用原公式即可运行：

mod_cox <- coxph(Surv(tstart, tstop, death==2) ~ sex + sex : rcs(alk.phos,4) + albumin + bili, 
                 id = factor(id), data = pbc2)

二、仅交互项模型的可视化 #

针对手动构造交互项的cph模型 #

利用rms::Predict直接指定alk.phos和sex的取值范围，工具会自动关联提前构造的变量计算预测值：

# 生成覆盖alk.phos全范围的预测数据，包含sex的两个水平
wadj <- rms::Predict(mod_adj, 
                     alk.phos = seq(min(pbc2$alk.phos, na.rm=TRUE), 
                                    max(pbc2$alk.phos, na.rm=TRUE), length.out=100), 
                     sex = c(0,1), 
                     conf.type='mean')
# 可视化
ggplot(wadj, aes(x=alk.phos, y=yhat, color=factor(sex), fill=factor(sex))) +
  geom_line(linewidth=1) +
  geom_ribbon(aes(ymin=lower, ymax=upper), alpha=0.2, color=NA) +
  labs(x="Alkaline Phosphatase", y="Predicted Hazard Ratio", color="Sex", fill="Sex") +
  theme_bw()

针对coxph模型 #

通过predict手动计算风险比及置信区间，再用ggplot绘制：

# 生成预测数据集：固定其他协变量为均值，覆盖alk.phos范围与sex水平
pred_data <- expand.grid(
  alk.phos = seq(min(pbc2$alk.phos, na.rm=TRUE), max(pbc2$alk.phos, na.rm=TRUE), length.out=100),
  sex = c(0,1),
  albumin = mean(pbc2$albumin, na.rm=TRUE),
  bili = mean(pbc2$bili, na.rm=TRUE),
  tstart = 0,
  tstop = max(pbc2$tstop, na.rm=TRUE)
)

# 计算风险比与标准误
pred_fit <- predict(mod_cox, newdata=pred_data, type="risk", se.fit=TRUE)
pred_data$risk <- pred_fit$fit
pred_data$se_risk <- pred_fit$se.fit

# 计算95%置信区间
pred_data$lower <- pred_data$risk * exp(-1.96 * pred_data$se_risk)
pred_data$upper <- pred_data$risk * exp(1.96 * pred_data$se_risk)

# 可视化
ggplot(pred_data, aes(x=alk.phos, y=risk, color=factor(sex), fill=factor(sex))) +
  geom_line(linewidth=1) +
  geom_ribbon(aes(ymin=lower, ymax=upper), alpha=0.2, color=NA) +
  labs(x="Alkaline Phosphatase", y="Predicted Hazard Ratio", color="Sex", fill="Sex") +
  theme_bw()

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