Loop in Cox regression

Question

I am trying to run a cox regression for 1000 variables (exposure) as below

varlist <- names(dataset)[275:1275]


sumtables <- lapply(varlist, function(i) {
     
      iformula <- as.formula(sprintf("Surv(time_cox, events) ~ %s + age +age2 ", i))      
      x <- coxph(iformula, data=dataset, na.action=na.omit)

      summary(x)[7][[1]] #####  summary(x)[8][[1]]
     
})

it works well, but I don't know how to extract the data (for each variable (beta and se)) and run the benjamini-hochberg on p-values. any help is appreciated! Thanks

Answer 1

You could use purrr::map to get a tidy dataframe of all your coefficients, se's and p values etc. from the vector of tested exposures. Modifying a little from your code above to work with veteran dataset:

library(survival)
library(tidyverse)

exp_vars <- names(veteran[, c(1, 2, 5, 6, 8)])

tibble(exp_vars) %>%
  group_by(exp_vars) %>%
  mutate(cox_mod = map(exp_vars, function(exposure) {
    iformula <-
      as.formula(sprintf("Surv(time, status) ~ %s + age", exposure))
    x <- coxph(iformula, data = veteran, na.action = na.omit)
    
    x
  }),
  coefs = list(rownames_to_column(data.frame(
    summary(cox_mod[[1]])$coefficients
  )))) %>%
  unnest(coefs)

#> # A tibble: 12 x 8
#> # Groups:   exp_vars [5]
#>    exp_vars cox_mod rowname             coef exp.coef. se.coef.       z Pr...z..
#>    <chr>    <list>  <chr>              <dbl>     <dbl>    <dbl>   <dbl>    <dbl>
#>  1 trt      <coxph> trt             -0.00365     0.996  0.183   -0.0200 9.84e- 1
#>  2 trt      <coxph> age              0.00753     1.01   0.00966  0.779  4.36e- 1
#>  3 celltype <coxph> celltypesmallc~  0.992       2.70   0.254    3.91   9.40e- 5
#>  4 celltype <coxph> celltypeadeno    1.16        3.17   0.293    3.94   8.07e- 5
#>  5 celltype <coxph> celltypelarge    0.235       1.27   0.278    0.848  3.97e- 1
#>  6 celltype <coxph> age              0.00590     1.01   0.00935  0.631  5.28e- 1
#>  7 karno    <coxph> karno           -0.0337      0.967  0.00520 -6.48   8.94e-11
#>  8 karno    <coxph> age             -0.00239     0.998  0.00908 -0.263  7.92e- 1
#>  9 diagtime <coxph> diagtime         0.00943     1.01   0.00892  1.06   2.90e- 1
#> 10 diagtime <coxph> age              0.00797     1.01   0.00961  0.830  4.07e- 1
#> 11 prior    <coxph> prior           -0.0135      0.987  0.0201  -0.674  5.00e- 1
#> 12 prior    <coxph> age              0.00715     1.01   0.00955  0.749  4.54e- 1

^{Created on 2022-03-16 by the reprex package (v2.0.1)}

Answer 2

I am assuming here that all the variables in varlist are either binary or numeric.

sumtables <- lapply(varlist, function(i) {

  iformula <- as.formula(sprintf("Surv(time_cox, events) ~ %s + age +age2 ", i))
  x <- coxph(iformula, data=dataset, na.action=na.omit)
  data.frame(pvalue = drop1(x, scope = i, test = "Chisq")[2,4],
             coef = coef(x)[i])
})