问题描述
由于我的真实数据违反了比例风险假设,我正在使用 AFT 模型,尝试计算感兴趣的研究组的调整后生存概率。下面的例子是关于 kidney 数据的,我试图遵循 ciTools vignette。
library(tidyverse)
library(ciTools)
library(here)
library(survival)
library(survminer)
#data
kidney
模型
fit1 = survreg(Surv(time,censored) ~ sex + age + disease,data = kidney)
Call:
survreg(formula = Surv(time,data = kidney)
Coefficients:
(Intercept) sexfemale age diseaseGN diseaseAN diseasePKD
8.41830937 -0.93959839 -0.01578812 -0.25274448 -0.38306425 -0.32830433
Scale= 1.642239
Loglik(model)= -122.1 Loglik(intercept only)= -122.7
Chisq= 1.33 on 5 degrees of freedom,p= 0.931
n= 76
添加至少存活 365 天的两性存活概率
probs = ciTools::add_probs(kidney,fit1,q = 365,name = c("prob","lcb","ucb"),comparison = ">")
probs
试图绘制两性的一年生存概率,但 geom_point 有多个点估计?
在我看来,这些点估计值是针对每个 年龄 值给出的。我可以编辑预测,使其针对平均或中位年龄吗?
probs %>% ggplot(aes(x = prob,y = sex)) +
ggtitle("1-year survival probability") +
xlim(c(0,1)) +
theme_bw() +
geom_point(aes(x = prob),alpha = 0.5,colour = "red")+
geom_linerange(aes(xmin = lcb,xmax = ucb),alpha = 0.5)
fit2 = survreg(Surv(time,censored) ~ sex,data = kidney)
probs2 = ciTools::add_probs(kidney,fit2,comparison = ">")
probs2 %>% ggplot(aes(x = prob,alpha = 0.5)
问题:
- 如何获得两性的调整后生存概率?或者,如果这是不可能的,那么可能的替代方案是什么?代码将有助于替代。
- 如果我想获得两性和不同时间点的调整后生存概率,我应该在 ciTools::add_probs() 函数中编辑“q”值吗?例如:q = 30 一个月; q = 90 三个月等。或者我应该为每个时间段运行一个单独的模型?
解决方法
根据您设置这些模型的方式,会根据模型中包含的协变量值为 kidney
数据集中的每个人返回预测。
在您的第一个模型中,您已包含 sex + age + disease
,以便您对数据集中这 3 个协变量值的每个组合进行预测。
在第二个模型中,您只包含 sex
作为预测变量,因此您只能获得基于性别的预测。
生存模型预测函数允许您指定一组协变量值,以便在新数据框中进行预测。根据 add_probs.survreg
的 manual,您可以通过在函数的 df
参数中指定具有指定协变量值的新数据框来实现。您在那里使用了 kidney
数据框,因此您可以预测所有这些情况。
我对 ciTools::add_probs
不太熟悉,但此类软件通常会(没有警告)接受您指定的协变量的值,然后对您不知道的协变量使用某种类型的“平均”值' t 指定。由于“平均值”对于 disease
等分类协变量没有太大意义,因此通常最好为所有协变量指定一组完整的有用值。
rms package
in R 中的函数通常在提供有用的预测方面做得更好,因为它们根据 {{1 }} 函数,然后将其输出指定为系统选项。此软件包的学习曲线有点陡峭,但如果您要进行大量生存或其他回归建模,则非常值得。