5

我有以下模型：

coxph(Surv(fulength, mortality == 1) ~ pspline(predictor))

其中 fulength 是随访的持续时间（包括死亡率），predictor 是死亡率的预测因子。

上面命令的输出是这样的：

                         coef  se(coef) se2    Chisq DF   p    
pspline(predictor), line 0.174 0.0563   0.0562 9.52  1.00 0.002
pspline(predictor), nonl                       4.74  3.09 0.200

如何绘制此模型，以便在 y 轴上获得具有 95% 置信带和风险比的漂亮曲线？我的目标与此类似：

在此处输入图像描述

4

1 回答 1

5

这是当您在 rms-package 的 ?cph 中运行第一个示例时得到的：

n <- 1000
set.seed(731)
age <- 50 + 12*rnorm(n)
label(age) <- "Age"
sex <- factor(sample(c('Male','Female'), n, 
              rep=TRUE, prob=c(.6, .4)))
cens <- 15*runif(n)
h <- .02*exp(.04*(age-50)+.8*(sex=='Female'))
dt <- -log(runif(n))/h
label(dt) <- 'Follow-up Time'
e <- ifelse(dt <= cens,1,0)
dt <- pmin(dt, cens)
units(dt) <- "Year"
dd <- datadist(age, sex)
options(datadist='dd')
S <- Surv(dt,e)

f <- cph(S ~ rcs(age,4) + sex, x=TRUE, y=TRUE)
cox.zph(f, "rank")             # tests of PH
anova(f)
plot(Predict(f, age, sex)) # plot age effect, 2 curves for 2 sexes

在此处输入图像描述

因为 rms/Hmisc 包组合使用格点图，所以需要使用格点函数来完成具有边际年龄密度特征的注释。另一方面，如果您想将响应值更改为相对危险，您只需在 Predict 调用中添加一个 'fun=exp' 参数并关联图表即可获得：

png(); plot(Predict(f, age, sex, fun=exp), ylab="Relative Hazard");dev.off()

在此处输入图像描述

于 2015-02-07T19:48:19.193 回答