如何在 y 轴上绘制带有危险无线电、死亡概率或死亡率的受限三次样条曲线？答案

【问题标题】：How to plot restricted cubic spline with hazard radio, probability of mortality, or mortality rate on y-axis?如何在 y 轴上绘制带有危险无线电、死亡概率或死亡率的受限三次样条曲线？
【发布时间】：2015-07-21 22:44:34
【问题描述】：

1) 如何在以下示例中将 y 轴更改为“适合”的“优势比”、“死亡率概率”和“死亡率”？

2) 在以下示例中，如何将“fit2”的 y 轴更改为“风险比”？

library(Hmisc)
library(survival)
library(rms)

data(pbc)
d <- pbc
rm(pbc)
d$died <- ifelse(d$status == 2, 1, 0)
d$status <- ifelse(d$status != 0, 1, 0)

ddist <- datadist(d)
options(datadist='ddist')

fit <- lrm(status ~ rcs(age, 4), data=d)
(an <- anova(fit))
plot(Predict(fit), anova=an, pval=TRUE)

fit2 <- cph(Surv(time, status) ~  rcs(age, 4), data=d)
(an2 <- anova(fit2))
plot(Predict(fit2), anova=an, pval=TRUE)

期待您的帮助！

更新 1 按照 BondedDust 的回答，我写了以下内容：

# probability
getProbability <- function(x) {
     exp(x)/(1+exp(x))*100
}

fit <- lrm(status ~ rcs(age, 4), data=d)
(an <- anova(fit))
plot(Predict(fit, fun=getProbability), anova=an, pval=TRUE, ylab="Probability of death [%]")

# overall probability to die
table(d$status)
round(table(d$status)[[2]]/sum(table(d$status))*100, digits=1) # = 44.5%

由于死亡的总体概率是 44.5%，所以作为非统计学家的我，预测概率的计算和结果图似乎是正确的，不是吗？

【问题讨论】：

标签： r plot logistic-regression spline cox-regression

【解决方案1】：

如果你想要优势比，那么你需要添加一个fun=-argument 来转换为优势比比例：

plot(Predict(fit,fun=exp), anova=an, pval=TRUE, ylab="Odds ratio")

我不确定我知道你所说的changing to the "probability of mortality", and "mortality rate" for "fit" 是什么意思。逆 logit 函数是 exp(x)/(1+exp(x)) 但为了从系数构建事件或速率的估计值，您需要合并截距项。也许您可以从 fit 对象的拟合值组件中提取一些有用的东西来满足您的作业问题的要求。

> names(fit)
 [1] "freq"              "sumwty"            "stats"             "fail"             
 [5] "coefficients"      "var"               "u"                 "deviance"         
 [9] "est"               "non.slopes"        "linear.predictors" "penalty.matrix"   
[13] "info.matrix"       "weights"           "call"              "Design"           
[17] "scale.pred"        "terms"             "assign"            "na.action"        
[21] "fail"              "interceptRef"      "nstrata"          
> str(fit$linear.predictors)
 Named num [1:418] -0.187 -0.235 0.41 -0.24 -0.538 ...
 - attr(*, "names")= chr [1:418] "1" "2" "3" "4" ...

与转换对数赔率系数与赔率比系数相同的方法适用于将对数风险转换为风险比：

plot(Predict(fit2, fun=exp), anova=an, pval=TRUE, ylab="Hazard ratio")

【讨论】：