Fine Gray模型非线性/交互项问题：logHR绘图与10年风险预测

你当前用Fine Gray模型建模癌症特异性死亡风险（结局1），竞争风险为其他原因死亡（结局2），纳入了TNM分期、Age_wins的限制性立方样条（RCS）以及Age_wins与Disease的交互项。针对你的两个问题，解决方案如下：

一、绘制Age_wins与癌症死亡的logHR非线性关联（分亚型） #

cmprsk包没有类似rms包Predict()的一键函数，需要手动构造预测数据并计算logHR及置信区间：

步骤与代码 #

# 加载所需包
library(cmprsk)
library(rms)
library(ggplot2)
library(dplyr)

# 生成示例数据（复用你的代码）
set.seed(100)
data1 <- data.frame(Age_wins = runif(1000, min=18, max=80), 
                    Disease=sample(c(0,1), 1000, replace=TRUE, prob=c(0.7, 0.3)), 
                    TNM=factor(sample(1:4, 1000, replace=TRUE)), 
                    Outcome=sample(c(0, 1, 2), 1000, replace=TRUE, prob=c(0.7, 0.1, 0.2)), 
                    Followup=runif(1000, min=1, max=240)
)

# 拟合Fine Gray模型（避免attach，直接引用数据框）
cov1 <- model.matrix(~rcs(Age_wins,3)*Disease+TNM, data=data1)[,-1]
crr1 <- crr(ftime=data1$Followup, fstatus=data1$Outcome, cov1=cov1, failcode=1, cencode=0)

# 1. 构造预测数据集：覆盖全年龄范围、两种亚型，固定TNM为参考水平（TNM=1）
pred_data <- expand.grid(
  Age_wins = seq(18, 80, length.out = 100),  # 生成100个年龄点保证曲线平滑
  Disease = c(0, 1),
  TNM = factor(1, levels = 1:4)
)

# 2. 转换为模型要求的矩阵格式（与原模型结构完全一致）
pred_cov <- model.matrix(~rcs(Age_wins,3)*Disease+TNM, data=pred_data)[,-1]

# 3. 预测logHR（线性预测值）及标准误
pred_lp <- predict(crr1, newdata=pred_cov, type="lp", se.fit=TRUE)

# 4. 合并结果并计算95%置信区间
pred_data <- pred_data %>%
  mutate(
    logHR = pred_lp$fit,
    se = pred_lp$se.fit,
    lower = logHR - 1.96 * se,
    upper = logHR + 1.96 * se
  )

# 5. 绘制曲线
ggplot(pred_data, aes(x=Age_wins, y=logHR, color=factor(Disease))) +
  geom_line(linewidth=1) +
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=factor(Disease)), alpha=0.2, color=NA) +
  labs(x="诊断年龄", y="logHR（相对于参考组）", color="癌症亚型", fill="癌症亚型") +
  theme_bw() +
  geom_hline(yintercept=0, linetype="dashed", color="gray50")

说明 #

• 参考组为：TNM=1、Disease=0、Age_wins取RCS默认参考点（通常是变量中位数）
• type="lp"返回线性预测值，即logHR（相对于参考组的风险对数比）
• 置信区间用标准误×1.96计算，对应95%置信水平

二、计算不同亚型、TNM分期的10年绝对风险 #

cmprsk的predict()函数可以直接指定时间点（120个月），返回该时间点的累积竞争风险（即癌症特异性死亡的绝对风险），无需手动筛选数据行。

步骤与代码 #

# 1. 构造全组合预测数据集：覆盖全年龄、两种亚型、所有TNM分期
risk_pred_data <- expand.grid(
  Age_wins = seq(18, 80, length.out = 50),  # 可调整点数平衡平滑度与计算速度
  Disease = c(0, 1),
  TNM = factor(1:4, levels = 1:4)
)

# 2. 转换为模型矩阵
risk_pred_cov <- model.matrix(~rcs(Age_wins,3)*Disease+TNM, data=risk_pred_data)[,-1]

# 3. 预测10年（120个月）的累积绝对风险
risk_pred <- predict(crr1, newdata=risk_pred_cov, type="cumhaz", ftime=120)

# 4. 合并结果
risk_pred_data <- risk_pred_data %>%
  mutate(
    `10年癌症死亡绝对风险` = risk_pred
  )

# 查看结果示例
head(risk_pred_data)

# 可选：绘制不同TNM分期、亚型的风险曲线
ggplot(risk_pred_data, aes(x=Age_wins, y=`10年癌症死亡绝对风险`, color=factor(Disease), linetype=TNM)) +
  geom_line(linewidth=1) +
  labs(x="诊断年龄", y="10年癌症死亡绝对风险", color="癌症亚型", linetype="TNM分期") +
  theme_bw()

说明 #

• type="cumhaz"返回Fine Gray模型的累积风险函数值，即考虑竞争风险下的癌症特异性死亡概率
• ftime=120直接指定预测的时间点为120个月（10年），无需筛选原始数据中的随访时间
• 结果包含所有你需要的组合：两种亚型、4种TNM分期、全年龄范围的10年绝对风险

额外注意事项 #

• 避免使用attach(data1)，容易引发变量命名冲突，建议直接通过数据框引用变量
• RCS节点数设置为3是基础配置，若样本量足够，可尝试增加到4或5以捕捉更复杂的非线性关系
• 预测时必须保证模型矩阵与原模型结构完全一致，使用同一公式生成model.matrix是最稳妥的方式

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