这是当您运行 rms 包的 ?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()