我认为你在 rms 包的 rcs 上走在了正确的轨道上。其实rms自带了自己的Coxph版本，叫做cph。

您可能希望尝试以下方法

fit <- cph (Surv(survival, event) ~ statin + rcs(bloodpressure, 3), x=TRUE, y=TRUE, data=data)
# x, y for predict, validate and calibrate;
plot(Predict(fit), data=data)

您可以在 Harrell 教授（RMS 软件包的作者和同名书籍的作者）的讲稿中了解更多相关信息http://biostat.mc.vanderbilt.edu/wiki/pub/Main/RmS/rms.pdf至第 18.4 节

然而，这是需要考虑的另一面 - 您是否在模型中包含了所有相关变量？模型中不必只有一个变量来可视化其效果。指定不足的模型会产生有偏估计，因为它不能控制所有潜在的混杂变量。（参见例如：http ://en.wikipedia.org/wiki/Omitted-variable_bias ）

已编辑：要使 plot.Predict 方法正常工作，您需要在 cph 行之前添加以下两行。（我不得不承认我不知道这个的确切含义，但我在 Harrell 的笔记 18.1 中看到了这一点，它有助于为我解决错误消息）希望这有帮助~

dd <- datadist(data)
options(datadist = 'dd')

不要接受这个答案，它只是为了多样化：

在肿瘤学中，有一种更接近比例风险的不同方法。它的工作原理如下：将血压量表分成小间隔，并“局部”估计风险比。您将获得血压与风险比的关系图。Bonetti 和 Gelber (2004) 将其称为 STEPP。还有一个 R 包。但是您必须记住这种方法的缺点：它需要相当大的样本量，并且结果取决于间隔的（任意）长度。毕竟，解释性分析比验证性分析更有用。此外，您不会获得最佳血压，而只会获得一个间隔。唯一的优点是您不必指定危险函数的形状。

（这也是一个缺点，因为它会诱使您在分析之前少思考）

我会尝试包cox.ph中的族mgcv，它可用于gam函数实现的广义加法模型。广义加性模型 (GAM) 是样条模型的多元线性回归模拟。如果你有一个大数据集，不用担心，因为mgcv它也很快。拟合灵活的半参数 GAM 后，您可以观察曲线，然后开始考虑一个能很好地捕捉其形状的参数广义非线性模型。

对于 BP 和事件之间的关系的图形显示，您可以使用 ggplot2 和 R。这是使用您自己的示例数据的示例：

library(ggplot2)
ggplot(data, aes(x=bloodpressure, y=event))+stat_smooth()

阴影区域表示 95% 的置信区间。

根据他汀类药物的使用，可以获得两条曲线：

ggplot(data, aes(x=bloodpressure, y=event, group=factor(statin), color=factor(statin)))+stat_smooth()

这里服用他汀类药物的患者很少，因此曲线不完整。

上面的代码使用 'loess' 作为方法。由于您期望它是一条曲线，因此可以将代码更改如下：

ggplot(data, aes(x=bloodpressure, y=event))+stat_smooth(method = "lm", formula = y ~ poly(x, 2), size = 1)

您还可以显示血压组中的实际事件发生率。假设您将所有人分为 10 个血压组。对于大样本量，可以增加组数以获得更准确的曲线。以下曲线中的 y 轴给出了该组中发生事件的患者比例。

data$grp = cut(data$bloodpressure, 10)
aa = aggregate(event~grp, data=data, mean)
ggplot(aa, aes(x=grp, y=event))+geom_line(aes(group=1))

可以添加错误栏：

ddt = data.table(data)
aaa = ddt[,list(meanevent=mean(event), se_event=se(event)),by=grp]
ggplot(aaa, aes(x=grp, y=meanevent))+geom_line(aes(group=1)) + geom_errorbar(aes(ymin = meanevent-se_event, ymax = meanevent+se_event), width = 0.2)