r - 如何从 R 中的 coxph 和脆弱模型更改风险比的单位增量？

Question

我运行了一个coxph模型和一个脆弱模型，但现在我想将连续变量（年龄）的风险比更改为以 5 个单位的增量而不是 1 个单位的形式显示。R中是否有可以执行此类任务的函数？如果是这样，该功能是否也适用于脆弱模式？我用了这个包frailtypack。

library('survival')
data(veteran)

cox <- coxph(Surv(time, status) ~ age, data = veteran)
summary(cox)

# Call:
#   coxph(formula = Surv(time, status) ~ age, data = veteran)
# 
#   n= 137, number of events= 128 
# 
#         coef exp(coef) se(coef)     z Pr(>|z|)
# age 0.007500  1.007528 0.009565 0.784    0.433
# 
#     exp(coef) exp(-coef) lower .95 upper .95
# age     1.008     0.9925    0.9888     1.027
# 
# Concordance= 0.515  (se = 0.029 )
# Likelihood ratio test= 0.63  on 1 df,   p=0.4
# Wald test            = 0.61  on 1 df,   p=0.4
# Score (logrank) test = 0.62  on 1 df,   p=0.4

标签： r

Answer 1

@milan 的帖子回答了一个类似的问题，但不是所问的问题。由于年龄被分成几十年并建模为连续变量，因此风险比将比较受试者的年龄十年与下一个最年轻的十年。也就是说，51 岁与 49 岁或 59 岁与 41 岁的受试者的 HR 将相同，尽管他们之间相隔 2 年或 18 年。

无论如何，您建议的默认值是连续变量中的 1 个单位增量，在这种情况下是年龄。通过 1 个单位的变化比较主题并不总是有用的，尤其是当范围变得更大时。

您可以执行以下对模型来说很幼稚的操作，因此这应该适用于 a lm, glm, survival::coxph,frailtypack::frailtyPenal等。

library('survival')
data(veteran)

## 1-year increase in age
cox <- coxph(Surv(time, status) ~ age, data = veteran)
exp(coef(cox))
#      age
# 1.007528

对于像 Cox 回归这样的乘法模型，您可以在模型拟合后获得 x 单位的变化：

## 5-year increase in age
exp(coef(cox)) ^ 5
#      age
# 1.038211

## or equivalently
exp(coef(cox) * 5)
#      age
# 1.038211

但是，为年龄转换创建变量然后拟合模型更容易：

## or you can create a variable to model
veteran <- within(veteran, {
  age5 <- age / 5
})
cox5_1 <- coxph(Surv(time, status) ~ age5, data = veteran)
exp(coef(cox5_1))
#    age10
# 1.038211

cox5_2 <- coxph(Surv(time, status) ~ I(age / 5), data = veteran)
exp(coef(cox5_2))
# I(age/5)
# 1.038211

注意这里需要I在公式界面中使用，因为有些运算符在公式中有特殊含义。例如，lm(mpg ~ wt - 1, mtcars)和lm(mpg ~ I(wt - 1), mtcars)是两个不同的模型。

您可以在其他模型中使用这些方法，例如，frailtyPenal如果这确实是您正在使用的模型：

library('frailtypack')
fp <- frailtyPenal(Surv(time, status) ~ age, data = veteran, n.knots = 12, kappa = 1e5)
exp(fp$coef)
exp(fp$coef) ^ 5

fp5_1 <- frailtyPenal(Surv(time, status) ~ age5, data = veteran, n.knots = 12, kappa = 1e5)
fp5_2 <- frailtyPenal(Surv(time, status) ~ I(age / 5), data = veteran, n.knots = 12, kappa = 1e5)
exp(fp5_1$coef)
exp(fp5_2$coef)