Search code examples
tidyversecox-regressionbroom

R: tidy_add_n() doesn't work when looping over many columns for Cox models with reformulate/as.formula


I would like to run a loop for Cox models that adjust for different subsets of covariates. The desired output should include number of events, person-years, hazards ratios, and 95%CI. If I call coxph directly, tidy_add_n works perfectly, but when I use either reformulate or as.formula in the loop, the function does not create the two columns n_event and exposure but a new column named n which contains only missing value. Any ideas or suggestions to fix this issue or to work around would be very much appreciated. Thank you!

Below are the reproducible codes.

test1 <- list(time = c(4,3,1,1,2,2,3), 
              status = c(1,1,1,0,1,1,0), 
              x1 = as.factor(c(0,2,1,1,1,0,0)), 
              x2 = as.factor(c(0,0,0,0,1,1,1)))

#the function works if I call coxph directly 
coxph(Surv(time, status) ~ x1 + x2, test1) %>%
tidy_and_attach()%>%
tidy_add_reference_rows() %>%
tidy_add_estimate_to_reference_rows() %>%
tidy_add_header_rows() %>%
tidy_add_n()


#tidy_add_n does not work when I use reformulate or as.formula.
surv <- "Surv(time, status)"
c("x1","x2")%>%
reformulate(response = surv)%>%
coxph(data = test1)%>%
tidy_and_attach()%>%
tidy_add_reference_rows() %>%
tidy_add_estimate_to_reference_rows() %>%
tidy_add_header_rows() %>%
tidy_add_n()

Solution

  • I don't know why, but this has something to do with the magrittr pipe %>% vs the new native pipe |>.

    model <- coxph(Surv(time, status) ~ x1 + x2, test1) %>%
      tidy_and_attach()%>%
      tidy_add_reference_rows() %>%
      tidy_add_estimate_to_reference_rows() %>%
      tidy_add_header_rows() %>%
      tidy_add_n()
    
    
    surv <- "Surv(time, status)"
    model2 <- c("x1","x2") |> 
      reformulate(response = surv)  |> 
      coxph(data = test1) |> 
      tidy_and_attach() |> 
      tidy_add_reference_rows() |> 
      tidy_add_estimate_to_reference_rows()  |> 
      tidy_add_header_rows() |> 
      tidy_add_n()
    
    all.equal(model %>% data.frame(), model2 %>% data.frame())
    

    You get the output:

    [1] TRUE
    

    (Converting both to data.frame() removes the attributes, which are still different, for some reason, despite the resulting tibbles looking identical in the console).

    Someone with more R knowledge might find this helpful in diagnosing why this is happening. But if you just want a fix - use the new pipe.

    EDIT:

    It's only the output of reformulate that needs changing to the native pipe for the "fix" to work:

    model2 <- c("x1","x2") %>% 
      reformulate(response = surv) |>  
      coxph(data = test1) %>%
      tidy_and_attach() %>%
      tidy_add_reference_rows() %>%
      tidy_add_estimate_to_reference_rows()  %>% 
      tidy_add_header_rows() %>%
      tidy_add_n()
    
    
    all.equal(model %>% data.frame(), model2 %>% data.frame())
    
    [1] TRUE