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

Question

I ran some code below that looks at running Cox regression across multiple outcome types (stroke, cancer, respiratory) that appear in separate columns. purrr seems to do this quite well. But I would also like to

1. print the name of each outcome type above the corresponding regression model and
2. print the coefficients as hazard ratios with 95% CIs.

I know this is quite a big ask but is important since my real dataset has almost 20 outcome types. Any help would be much appreciated!

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))

Answer 1

I'm not quite sure if this is what you are looking for, but if you run the following code:

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)
})

It will print out:

#> 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)

And your list of models will be stored with the correct call above them:

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