在 coxph 输出和求幂系数上方打印结果名称(R 版本 4.1.2 (2021-11-01) -- "Bird Hippie")

Printing name of outcome above coxph output and exponentiating coefficients (R version 4.1.2 (2021-11-01) -- "Bird Hippie")

我 运行 下面的一些代码着眼于 运行 出现在单独列中的多种结果类型(中风、癌症、呼吸系统)的 Cox 回归。 purrr 似乎做得很好。不过我也想

  1. 在相应的回归模型上方打印每个结果类型的名称,并且
  2. 将系数打印为具有 95% 置信区间的风险比。

我知道这是一个很大的问题,但很重要,因为我的真实数据集有近 20 种结果类型。任何帮助将不胜感激!

library(survival)
library(purrr)

mydata <- read.table(header=T, 
text="age    Sex    survival    stroke cancer respiratory
51   2   1.419178082 2 1 1
60    1   5   1 2 2
49    2   1.082191781 2 2 2
83    1   0.038356164 1 1 2
68    2   0.77260274  2 1 2
44    2   2.336986301 1 2 1
76    1   1.271232877 1 2 2")

outcomes <- names(mydata[4:6])

purrr::map(outcomes, ~coxph(as.formula(paste("Surv(survival,", .x, ") ~  Sex + age")), 
mydata))

我不太确定这是否是您要查找的内容,但如果您运行以下代码:

result <- purrr::map(outcomes, function(x) {
  f <- as.formula(paste("Surv(survival,", x, ") ~  Sex + age"))
  model <- coxph(f, mydata)
  model$call$formula <- f
  s <- summary(model)
  cat(x, ':\n', paste0(apply(s$coefficients, 1, 
          function(x) {
            paste0("HR : ", round(exp(x[1]), 2), 
                   ' (95% CI ', round(exp(x[1] - 1.96 * x[3]), 2),
                   ' - ', round(exp(x[1] + 1.96 * x[3]), 2), ')')}),
        collapse = '\n'), '\n\n', sep = '')
  invisible(model)
})

它会打印出:

#> stroke:
#> HR : 650273590159.06 (95% CI 0 - Inf)
#> HR : 1.36 (95% CI 0.75 - 2.49)
#>
#> cancer:
#> HR : 1121.58 (95% CI 0 - 770170911.09)
#> HR : 1.33 (95% CI 0.78 - 2.28)
#>
#> respiratory:
#> HR : 24.1 (95% CI 0.31 - 1884.85)
#> HR : 1.2 (95% CI 0.99 - 1.45)

您的模型列表将与上面的正确调用一起存储:

result
#> [[1]]
#> Call:
#> coxph(formula = Surv(survival, stroke) ~ Sex + age, data = mydata)
#> 
#>          coef exp(coef)  se(coef)     z     p
#> Sex 2.720e+01 6.503e+11 2.111e+04 0.001 0.999
#> age 3.105e-01 1.364e+00 3.066e-01 1.013 0.311
#> 
#> Likelihood ratio test=6.52  on 2 df, p=0.03834
#> n= 7, number of events= 3 
#> 
#> [[2]]
#> Call:
#> coxph(formula = Surv(survival, cancer) ~ Sex + age, data = mydata)
#> 
#>          coef exp(coef)  se(coef)     z     p
#> Sex    7.0225 1121.5843    6.8570 1.024 0.306
#> age    0.2870    1.3325    0.2739 1.048 0.295
#> 
#> Likelihood ratio test=2.58  on 2 df, p=0.2753
#> n= 7, number of events= 4 
#> 
#> [[3]]
#> Call:
#> coxph(formula = Surv(survival, respiratory) ~ Sex + age, data = mydata)
#> 
#>         coef exp(coef) se(coef)     z      p
#> Sex  3.18232  24.10259  2.22413 1.431 0.1525
#> age  0.18078   1.19815  0.09772 1.850 0.0643
#> 
#> Likelihood ratio test=5.78  on 2 df, p=0.05552
#> n= 7, number of events= 5